COMMENT ⊗   VALID 00023 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00004 00002	%folio 397 galley 1 (C) Addison-Wesley 1978	-
C00018 00003	%folio 403 galley 2 (C) Addison-Wesley 1978	-
C00034 00004	%folio 407 galley 3 (C) Addison-Wesley 1978	¬
C00052 00005	%folio 411 galley 4a (C) Addison-Wesley 1978	¬
C00056 00006	%folio 412 galley 4b (C) Addison-Wesley 1978	-
C00069 00007	%folio 415 galley 1 (C) Addison-Wesley 1978	-
C00080 00008	%folio 417 galley 2 (C) Addison-Wesley 1978	-
C00092 00009	%folio 419 galley 3 (C) Addison-Wesley 1978	-
C00100 00010	%folio 421 galley 4 WARNING: Much tape unreadable! (C) Addison-Wesley 1978	-
C00109 00011	%folio 423 galley 5 (C) Addison-Wesley 1978	-
C00119 00012	%folio 426 galley 6 (C) Addison-Wesley 1978	-
C00133 00013	%folio 429 galley 7 (C) Addison-Wesley 1978	-
C00147 00014	%folio 432 galley 8 (C) Addison-Wesley 1978	-
C00159 00015	%folio 434 galley 9 (C) Addison-Wesley 1978	-
C00171 00016	%folio 439 galley 10 WARNING: Some bad spots. (C) Addison-Wesley 1978	-
C00180 00017	%folio 441 galley 11 (C) Addison-Wesley 1978	-
C00194 00018	%folio 444 galley 12 (C) Addison-Wesley 1978	-
C00206 00019	%folio 448 galley 13 (C) Addison-Wesley 1978	-
C00215 00020	%folio 451 galley 14 (C) Addison-Wesley 1978	-
C00225 00021	%folio 460 galley 15 (C) Addison-Wesley 1978	-
C00233 00022	%folio 461 galley 16 (C) Addison-Wesley 1978	-
C00245 00023	
C00246 ENDMK
C⊗;
%folio 397 galley 1 (C) Addison-Wesley 1978	-
\sectionbegin{4.4. RADIX CONVERSION}
I{\:cF MEN} had invented arithmetic by counting with
their two fists or their eight fingers, instead of their ten
``digits,'' we would never have to worry about writing binary-decimal
conversion routines. (And we would perhaps never have learned
as much about number systems.) In this section, we shall discuss
the conversion of numbers from positional notation with one
radix into positional notation with another radix; this process
is, of course, most important on binary computers when converting
decimal input data into binary form, and converting binary answers
into decimal form.

\subsectionbegin{A. The four basic methods} Binary-decimal
conversion is one of the most machine-dependent operations of
all, since engineers keep inventing different ways to provide
for it in the computer hardware. Therefore we shall discuss
only the general principles involved, from which a programmer
can select the procedure most well suited to his machine.

 We shall assume that only nonnegative numbers enter
into the conversion, since the manipulation of signs is easily
accounted for.

 Let us assume that we are converting from radix
$b$ to radix $B$. (The methods can also be generalized to mixed-radix
notations, as shown in exercises 1 and 2.)
Most radix-conversion routines
are based on multiplication and division, using one of the following
four schemes:

\yyskip
1) Conversion of integers (radix point at the right).

\yyskip
\noalign$\bullet$ Method (1a) {\sl Division by $B$} (using
radix-$b$ arithmetic).\xskip Given an integer number $u$, we obtain
its radix-$B$ representation $(U↓M \ldotsm U↓1U↓0)↓B$ as follows:
$$\baselineskip14pt
\eqalign{U↓0 ⊗= u \mod B\cr
U↓1 ⊗= \lfloor u/B\rfloor \mod B\cr
U↓2 ⊗= \lfloor \lfloor u/B\rfloor /B\rfloor \mod B\cr
\ldots\cr}$$
etc., stopping when $\lfloor \ldotsm \lfloor \lfloor
u/B\rfloor /B\rfloor \ldotsm /B\rfloor =0$.

\yyskip
\noalign$\bullet$ Method (1b) {\sl Multiplication by $b$} (using
radix-$B$ arithmetic).\xskip If $u$ has the radix-$b$ representation
$(u↓m \ldotsm u↓1u↓0)↓b$, we can use radix-$B$ arithmetic to
evaluate the polynomial $u↓mb↑m + \cdots + u↓1b + u↓0=u$ in
the form
$$\biglp (\ldotsm (u↓mb + u↓{m-1})b+\cdotss)b+u↓1\bigrp b+u↓0.$$

\yyskip 2) Conversion of fractions (radix
point at the left).\xskip Note that it is often impossible to express
a terminating radix-$b$ fraction $(0.u↓{-1}u↓{-2} \ldotsm u↓{-m})↓b$
{\sl exactly} as a terminating radix-$B$ fraction $(0.U↓{-1}U↓{-2}
\ldotsm U↓{-M})↓B$. For example, the fraction ${1\over 10}$ has
the infinite binary representation $(0.0 0011 0011 0011 . . .)↓2$.
Therefore methods of rounding the result to $M$ places are sometimes
of interest.

\yyskip\noindent$\bullet$ Method (2a) {\sl Multiplication by $B$} (using
radix-$b$ arithmetic).\xskip Given a
fractional number $u$, we obtain successive digits of its radix-$B$
representation
$(.U↓{-1}U↓{-2}\ldotsm)↓B$ as follows:
$$\baselineskip14pt
\eqalign{U↓{-1} ⊗= \lfloor uB\rfloor \cr
U↓{-2} ⊗= \lfloor \{uB\} B\rfloor \cr
U↓{-3} ⊗= \lfloor \{\{uB\}B\} B\rfloor \cr
\ldots\cr}$$
where $\{x\}$ denotes $x \mod 1 = x - \lfloor
x\rfloor $. If it is desired to round the result to $M$ places,
the computation can be stopped after $U↓{-M}$ has been calculated,
and $U↓{-M}$ should be increased by unity if $\{\ldotsm\{\{nB\}B\}\ldotsm
B\}$ is greater than ${1\over 2}$. (Note, however,
that this may cause a number of carries to occur, which must
be incorporated into the answer using radix-$B$ arithmetic.
It would be simpler to add the constant ${1\over 2}B↑{-M}$
to the original number $u$ before the calculation begins, but
this may lead to a terribly incorrect answer when ${1\over 2}B↑{-M}$
cannot be represented exactly as a radix-$b$ number inside the
computer. Note further that it is possible for the answer to
round up to $(1.00 \ldotsm 0)↓B$, if $b↑m ≥ 2B↑M$.)

Exercise 3 shows how to extend this method so that $M$ is {\sl variable},
just large enough to represent the original number to a specified
accuracy; in this case the problem of carries does not occur.

\yyskip\noindent$\bullet$ Method (2b) {\sl Division by $b$} (using
radix-$B$ arithmetic).\xskip If $u$ has the radix-$b$ representation
$(0.u↓{-1}u↓{-2} \ldotsm u↓{-m})↓b$, we can use radix-$B$ arithmetic
to evaluate $u↓{-1}b↑{-1} + u↓{-2}b↑{-2} + \cdots + u↓{-m}b↑{-m}$
in the form
$$\biglp(\ldotsm (u↓{-m}/b + u↓{1-m})/b + \cdots
+ u↓{-2})/b + u↓{-1}\bigrp /b.$$
Care should be taken to control errors
that might occur due to truncation or rounding in the division
by $b$; these are often negligible, but not always.

\yyskip
To summarize, Methods (1a), (1b), (2a), and (2b)
give us two choices for a conversion process, depending on whether
our number is an integer or a fraction. And it is certainly
possible to convert between integers and fractions by multiplying
or dividing by an appropriate power of $b$ or $B$; therefore
there are at least four methods to choose from when trying to
do a conversion.

\subsectionbegin{B. Single-precision conversion}
To illustrate these four methods, let us assume that \MIX\ is
a binary computer, and suppose that we want to convert a binary
integer $u$ to a decimal integer, that is, $b = 2$ and $B =
10$. Method (1a) could be programmed as follows:
$$\vcenter{\mixthree{\!
⊗ENT1⊗0⊗Set $j ← 0$.\cr
⊗LDX⊗U\cr
⊗ENTA⊗0⊗Set $\rAX ← u$.\cr
1H⊗DIV⊗=10=⊗$(\rA, \rX) ← (\lfloor \rAX/10\rfloor , \rAX \mod 10)$.\cr
⊗STX⊗ANSWER,1⊗$U↓j ←\rX$.\cr
⊗INC1⊗1⊗$j ← j + 1$.\cr
⊗SRAX⊗5⊗$\rAX ← \rA$.\cr
⊗JXP⊗1B⊗Repeat until result is zero.\quad\blackslug\cr}}\eqno (1)$$
This requires $18M + 4$ cycles to obtain $M$ digits.

The above method uses division by 10; Method (2a)
uses {\sl multiplication} by 10, so it might be a little faster.
In order to use Method (2a), we must deal with fractions, and
this leads to an interesting situation. Let $w$ be the word
size of the computer, and assume that $u < 10↑n < w$. With a
single division we can find $q$ and $r$, where
$$wu=10↑nq + r,\qquad 0 ≤ r < 10↑n.\eqno (2)$$
Now if we apply Method (2a) to the fraction $(q
+ 1)/w$, we will obtain the digits of $u$ from left to right,
in $n$ steps, since
$$\left\lfloor 10↑n {q + 1\over w}\right\rfloor
 = \left\lfloor u + {10↑n - r\over w}\right\rfloor= u.\eqno (3)$$
(This idea is due to P. A. Samet, {\sl Software---Practice
and Experience \bf 1} (1971), 93--96.)

Here is the corresponding \MIX\ program:
$$\vcenter{\mixthree{\!
⊗JOV⊗OFLO⊗Ensure overflow is off.\cr
⊗LDA⊗U\cr
⊗LDX⊗=$10↑n$=⊗$\rAX ← wu + 10↑n$.\cr
⊗DIV⊗=$10↑n=⊗$\rA ← q + 1$, $\rX ← r$.\cr
⊗JOV⊗ERROR⊗Jump if $u ≥ 10↑n$.\cr
⊗ENT1⊗$n$-1⊗Set $j ← n-1$.\cr
2H⊗MUL⊗=10=⊗Now imagine radix point at left, $\rA = x$.\cr
⊗STA⊗ANSWER,1⊗Set $U↓j ← \lfloor 10x\rfloor$.\cr
⊗SLAX⊗5⊗$x ← \{10x\}$.\cr
⊗DEC1⊗1⊗$j ← j - 1$.\cr
⊗J1NN⊗2B⊗Repeat for $n > j ≥ 0$.\quad\blackslug\cr}}\eqno (4)$$
This slightly longer routine requires $16n
+ 19$ cycles, so it is a little faster than (1) if $n = M ≥
8$; when leading zeros are present, (1) will be faster.

Program (4) as it stands cannot be used to convert
integers $u ≥ 10↑m$ when $10↑m < w < 10↑{m+1}$, since we would
need to take $n = m + 1$. In this case we can obtain the leading
digit of $u$ by computing $\lfloor u/10↑m\rfloor $; then $u
\mod 10↑m$ can be converted as above with $n = m$.

The fact that the answer digits are obtained from
left to right may be an advantage in some applications (e.g.,
when typing out the answer one digit at a time). Thus we see
that a fractional method can be used for conversion of integers,
although the use of inexact division makes a little bit of numerical
analysis necessary.

A modification of Method (1a) can be used to avoid
division by 10, by replacing it with two multiplications. It
is worth mentioning this modification here, because radix conversion
is often done by small ``satellite'' computers that have no
division capability. If we let $x$ be an approximation to ${1\over
10}$, so that ${1\over 10}$ < x < ${1\over 10} + 1/w$, it
is easy to prove (see exercise 7) that $\lfloor ux\rfloor =
\lfloor u/10\rfloor$ or $\lfloor u/10\rfloor + 1$, so long as
$0 ≤ u < w$. Therefore, if we compute $u - 10\lfloor ux\rfloor
$, we will be able to determine the value of $\lfloor u/10\rfloor$:
$$\lfloor u/10\rfloor = \left\{\vcenter{\halign{$\lft{#},$\qquad
⊗if $u-10\lfloor ux\rfloor#\cr
\lfloor ux\rfloor⊗≥0;\cr
\lfloor ux\rfloor-1⊗<0.\cr}}\right.\eqno(5)$$
This procedure simultaneously determines $u\mod
10$. A \MIX\ program for conversion using this idea appears in
exercise 8; it requires about 33 cycles per digit.
%folio 403 galley 2 (C) Addison-Wesley 1978	-
The procedure represented by (5) can be used effectively even
if the computer has no built-in multiplication instruction,
since multiplication by 10 consists of two shifts and one addition
($10 = 2↑3 + 2$), and multiplication by ${1\over 10}$ can also
be done by combining a few shifting and adding operations judiciously
(see exercise 9).

Method (1b) could also be used to convert from binary to decimal,
but to do this we need to simulate doubling in a {\sl decimal}
number system. This idea is generally most suitable for incorporation
into computer hardware; however, it is possible to program the
doubling process for decimal numbers, using binary addition,
shifting, and extraction (``logical \.{AND}'' on each bit in the
register), as shown in Table 1.

\topinsert{\tablehead{Table 1}
\vskip 3pt plus 2pt minus 1pt
\ctrline{DOUBLING A BINARY-CODED DECIMAL NUMBER}
\vskip 3pt plus 2pt minus 1pt
\eightpoint
\halign{\rt{#} ⊗\lft{#}\quad⊗$\ctr{#}$⊗ $\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗
$\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗
$\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗$\ctr{#}$⊗\quad\rt{#}⊗$\null#\hfill$\cr
⊗\sl Operation⊗⊗⊗⊗⊗⊗⊗⊗\hskip 0pt minus 100pt \sl General form\hskip 0pt minus 100pt
⊗⊗⊗⊗⊗⊗\sl Example\hfill\cr
\noalign{\vskip 2pt}
1.⊗Given⊗⊗u↓1⊗u↓2⊗u↓3⊗u↓4⊗u↓5⊗u↓6⊗u↓7⊗u↓8⊗u↓9⊗
u↓{10}⊗u↓{11}⊗u↓{12}⊗0011 0110 1001⊗= 3\;6\;9\cr
⊗number\cr
2.⊗Add 3 to⊗⊗v↓1⊗v↓2⊗v↓3⊗v↓4⊗v↓5⊗v↓6⊗v↓7⊗v↓8⊗v↓9⊗v↓{10}⊗
v↓{11}⊗v↓{12}⊗0110 1001 1100\cr
⊗each digit\cr
3.⊗Shift left⊗v↓1⊗v↓2⊗v↓3⊗v↓4⊗v↓5⊗v↓6⊗v↓7⊗v↓8⊗v↓9⊗v↓{10}⊗v↓{11}⊗
v↓{12}⊗0⊗0 1101 0011 1000\cr
⊗one\cr
4.⊗Extract low⊗v↓1⊗0⊗0⊗0⊗v↓5⊗0⊗0⊗0⊗v↓9⊗0⊗0⊗0⊗0⊗0 0001 0001
0000\cr
⊗bit\cr
5.⊗Shift right⊗⊗0⊗v↓1⊗0⊗0⊗0⊗v↓5⊗0⊗0⊗0⊗v↓9⊗0⊗0⊗
0000 0100 0100\cr
⊗two\cr
6.⊗Shift right⊗⊗0⊗v↓1⊗v↓1⊗0⊗0⊗v↓5⊗v↓5⊗0⊗0⊗v↓9⊗v↓9⊗0⊗
0000 0110 0110\cr
⊗one and add\cr
7.⊗Add result⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗\ast⊗0
⊗0 1101 1001 1110\cr
⊗of step 3\cr
8.⊗Subtract 6⊗y↓1⊗y↓2⊗y↓3⊗y↓4⊗y↓5⊗y↓6⊗y↓7⊗y↓8⊗y↓9⊗y↓{10}⊗y↓{11}⊗
y↓{12}⊗0⊗0 0111 0011 1000⊗= 7\;3\;8\cr
⊗from each\cr}}

This method changes each individual digit $d$
into $\biglp (d + 3) \times 2 + 0) - 6 = 2d$ when $0 ≤ d ≤ 4$,
and into $\biglp (d + 3) \times 2 + 6\bigrp - 6 = (2d - 10)
+ 2↑4$ when $5 ≤ d ≤ q$; and that is just what is needed to
double decimal numbers encoded with 4 bits per digit.

Another relaed idea is to keep a table
of the powers of two in decimal form, and to add the appropriate
powers together by simulating decimal addition. A survey of
such bit-manipulation techniques appears in Section 7.1.

Finally, even Method (2b) can be used for the conversion of
binary integers to decimal integers. We can find $q$ as in
(2), and then we can simulate the decimal division of $q + 1$
by $w$, using a ``halving'' process (exercise 10) that is similar
to the doubling process just described, and retaining only the
first $n$ digits to the right of the radix point in the answer.
In this situation, Method (2b) does not seem to offer advantages
over the other three methods already discussed, but we have
confirmed the remark made earlier that at least four distinct
methods are available for converting integers from one radix
to another.

Now let us consider decimal-to-binary conversion ($b=10$, $B=2$).
Method (1a) simulates a decimal division by 2; this is feasible
(see exercise 10), but it is primarily suitable for incorporation
in hardware instead of programs.

Method (1b) is the most practical method
for decimal-to-binary conversion in the great majority of cases.
Here it is in \MIX\ code, assuming that there are at least two
digits in the number $(u↓m \ldotsm u↓1u↓0)↓{10}$ being converted:
$$\vcenter{\mixthree{\!
⊗ENT1⊗M-1⊗Set $j ← m - 1$.\cr
⊗LDA⊗INPUT+M⊗Set $U ← u↓m$.\cr
1H⊗MUL⊗=10=\cr
⊗SLAX⊗5\cr
⊗ADD⊗INPUT,1⊗$U ← 10U + u↓j$.\cr
⊗DEC1⊗1\cr
⊗J1NN⊗1B⊗Repeat for $m > j ≥ 0$.\quad\blackslug\cr}}\eqno (6)$$
Note again that adding and shifting may be
substituted for multiplication by 10.

A trickier but perhaps faster method, which uses about $\lg m$
multiplications, extractions, and additions instead of $m$ multiplications
and additions, is described in exercise 19.

For the conversion of decimal fractions $(0.u↓{-1}u↓{-2} \ldotsm
u↓{-m})↓{10}$, we can use Method (2b), or, more commonly, we
can convert the integer $(u↓{-1}u↓{-2} \ldotsm u↓{-m})↓{10}$
by Method (1b) and then divide by $10↑m$.

\subsectionbegin{C. Hand calculation} It is occasionally
necessary for computer programmers to convert numbers
by hand, and since this is a subject not yet taught in elementary
schools, it may be worth while to examine it briefly here. There
are very simple pencil-and-paper methods for converting between
decimal and octal notations, and these methods are easily learned,
so they ought to be more widely known.

\noindent {\sl Converting octal integers to decimal. \xskip}
The simplest conversion is from octal to decimal; this technique
was apparently first published by Walter Soden, {\sl Math.\ Comp.\
\bf 7} (1953), 273--274. To do the conversion, write down the
given octal number; then at the $k$th step, double the $k$ leading
digits using decimal arithmetic, and subtract this from the
$k + 1$ leading digits using decimal arithmetic. The process
terminates in $n - 1$ steps if the given number has $n$ digits.
It is a good idea to insert a radix point to show which digits
are being doubled, as shown in the following example, in order
to prevent embarrassing mistakes.

{\baselineskip0pt\lineskip0pt\def\\{\lower 2.5pt\vjust to 11pt{}}
\lower 1pt\vjust to 2.4pt{}⊗\vjust{\hrule width 18pt}⊗\vjust{\hrule width 18pt}⊗
\vjust{\hrule width 18pt}⊗\vjust{\hrule width 36pt}\cr
\thbegin Example 1. Convert $(5325121)↓8$ to decimal.
$$\vjust{\baselineskip0pt\lineskip0pt
\halign{\hjust to 100pt{\hfill#}⊗\hjust to 100pt{#\hfill}⊗\hjust to 148pt
{#\hfill}\cr
\def\dot{\hjust to 5pt{\hfill.\hfill}
\def\\{\lower 2.5pt\vjust to 12pt{}}
\def\¬{$-$ \\}
\def\bar{\lower 1pt\vjust to 2.4pt{}}
\\⊗5\dot3\92\95\91\92\91\cr
\¬⊗1\90\cr
\bar⊗\vjust{\hrule width 15pt}\cr
\\4\93\dot2\95\91\92\91\cr
\¬⊗\9\98\96\cr
\bar⊗\vjust{\hrule width 25pt}\cr
|qleft ⊗3 4 6 . 5 1 2 1\cr
\hfill - 6 9 2\cr
\-2skip \quad ------------
|qleft ⊗2 7 7 3 . 1 2 1\cr
\hfill - 5 5 4 6\cr
\-2skip \quad ------------
|qleft ⊗2 2 1 8 5 . 2 1\cr
\hfill - 4 4 3 7 0\cr
\-2skip \quad ---------------
|qleft ⊗1 7 7 4 8 2 . 1\cr
\hfill - 3 5 4 9 6 4\cr
\-2skip \quad ---------------------
|qleft ⊗1 4 1 9 8 5 7\cr
⊗\qquad \qquad \qquad \quad {\sl Answer{\sl :} (1419857)↓{10}.\cr
\9skip} \quad A reasonably good check on the computations may
be had by ``casting out mines'': The sum of the digits of the
decimal number must equal the alternating sum and difference
of the digits of the octal number, with the rightmost digit
of the latter given a plus sign, except for a multiple of nine.
In the above example, we have 1 + 4 + 1 + 9 + 8 + 5 + 7 = 35,
and 1 - 2 + 1 - 5 + 2 - 3 + 5 = -1, and the difference is 36
(a multiple of 9). If this test fails, it can be applied to
the $k + 1$ leading digits after the $k$th step, and the error
can be located using a ``binary search'' procedure; i.e., start
by checking the middle result, then use the same procedure on
the first or second half of the calculation, depending on whether
the middle result is incorrect or correct.

Of course, the ``casting-out-nines'' process is only about 89
percent reliable; \quad Of course, the ``casting-out-nines''
process is only about 89 percent reliable; an even better check
is to convert the answer back to octal by using an inverse method,
which we will now consider.

|qleft \12skip {\sl Converting decimal integers to octal. \xskip}
A similar procedure can be used for the opposite conversion:
Write down the given decimal number; at the $k$th step, double
the $k$ leading digits using {\sl octal} arithmetic, and {\sl
add} these to the $k + 1$ leading digits using {\sl octal} arithmetic.
The process terminates in $n - 1$ steps if the given number
has $n$ digits.

|qleft \6skip {\bf Example 2. \xskip} Convert (1419857)$↓{10}
to octal$.

|qleft \6skip \qquad \qquad \qquad \qquad |tab 1 .}24 1 9 8
5 7⊗⊗\quad 2\cr
\-2skip ⊗------\cr
⊗1 6 . 1 9 8 5 7\cr
⊗\quad 3 4\cr
\-2skip \qquad \qquad \qquad \qquad ---------
|qleft ⊗2 1 5 . 9 8 5 7\cr
⊗\quad 4 3 2\cr
\-2skip \qquad \qquad \qquad \qquad ------------
|qleft ⊗2 6 1 3 . 8 5 7\cr
⊗\quad 5 4 2 6\cr
\-2skip \qquad \qquad \qquad \qquad ------------
|qleft ⊗3 3 5 6 6 . 5 7\cr
⊗\quad 6 7 3 5 4\cr
\qquad \qquad \qquad \qquad ---------------
|qleft ⊗4 2 5 2 4 1 . 7\cr
⊗1 0 5 2 5 0 2\cr
\-2skip \qquad \qquad \qquad \qquad ------------------
|qleft ⊗5 3 2 5 1 2 1\cr
\6skip ??M13}(Note that the nonoctal digits 8 and 9 center into
this octal cmputation.) The answer can be checked as discussed
above. This method was published by Charles P. Rozier, {\sl
IEEE Trans.\ \bf CE--11} (1962), 708--709.

|qleft \60skip {\sl Answer{\sl :} (5325121)↓8.}

|qleft \9skip |newtype \quad The two procedures just given are
essentially Method (1b) of the general radix-conversion procedures.
Doubling and subtracting in decimal notation is like multiplying
by 10 - 2 = 8; doubling and adding in octal notation is like
multiplying by 8 + 2 = 10. There is a similar method for hexadecimal/decimal
conversions, but it is a little more difficult since it involves
multiplication by 6 instead of by 2.

To keep these two methods straingt in our minds, it is not hard
to remember that we must subtract to go fr|newtype 58320---Computer
Programming\quad (Knuth/A
%folio 407 galley 3 (C) Addison-Wesley 1978	¬
ddison-Wesley)\quad f. 407\quad ch. 4\quad g. 3(c)

|qleft \12skip {\sl Converting fractions. \xskip} No equally
fast method of converting fractions manually is known; the best
way seems to be Method (2a), with doubling and adding or subtracting
to simplify the multiplications by 10 or by 8. In this case,
we reverse the addition-subtraction criterion, adding when we
convert to decimal and subtracting when we convert to octal;
we also use the radix of the given input number, {\sl not} the
radix of the answer, in this computation (see Examples 3 and
4). The process is about twice as hard as the above method for
integers.

|qleft \12skip ??M11.6}{\bf Example 3. \xskip} Convert (.14159)$↓{10}
to octal$.

|qleft \6skip . 1 4 1 5 9
2 8 3 1 8-
|qleft \-2skip ---------------
|qleft 1 . 1 3 2 7 2
|qleft \quad 2 6 5 4 4-
|qleft \-2skip \quad ---------------
|qleft \quad 1 . 0 6 1 7 6
|qleft \qquad \quad 1 2 3 5 2---
|qleft \-2skip \qquad ---------------
|qleft \qquad 0 . 4 9 4 0 8
|qleft \qquad \qquad 9 8 8 1 6-
|qleft \-2skip \qquad \quad ---------------
|qleft \qquad \quad 3 . 9 5 2 4 6
|qleft \qquad \quad \quad 1 9 0 5 2 8-
|qleft \-2skip \qquad \qquad ---------------
|qleft \qquad \qquad 7 . 6 2 1 1 2
|qleft \qquad \qquad \quad 1 2 4 2 2 4-
|qleft \qquad \qquad \quad ---------------
|qleft \qquad \qquad \quad 4 . 9 6 8 9 6

|qleft \6skip {\sl Answer{\sl :} (.110374 . . .)↓8.}

|qleft \6skip {\bf Example 4. \xskip} Convert (.110374)$↓8 to
decimal$.

|qleft \6skip . 1 1 0 3 7 4
|qleft \quad 2 2 0 7 7 0+
|qleft \-2skip ---------------
|qleft 1 . 3 2 4 7 3 0
|qleft \qquad 6 5 1 6 6 0+
|qleft \-2skip \quad ------------------
|qleft \quad 4 . 1 2 1 1 6 0
|qleft \qquad \quad 2 4 2 3 4 0+
|qleft \-2skip \qquad ------------------
|qleft \qquad 1 . 4 5 4 1 4 0
|qleft \qquad \quad 1 1 3 0 3 0 0+
|qleft \-2skip \qquad \quad ---------------------
|qleft \qquad \quad 5 . 6 7 1 7 0 0
|qleft \qquad \qquad 1 5 6 3 6 0 0+
|qleft \-2skip \qquad \qquad ---------------------
|qleft \qquad \qquad 8 . 5 0 2 6 0 0
|qleft \qquad \qquad \quad 1 2 0 5 4 0 0+
|qleft \-2skip \qquad \qquad \quad ---------------------
|qleft \qquad \qquad \quad 6 . 2 3 3 4 0 0

|qleft \6skip {\sl Answer{\sl :} (.141586 . . .)↓{10}.}

\subsectionbegin{D. Floating-point
conversion} When floating-point values are to be converted,
it is necessary to deal with both the exponent and the fraction
parts simultaneously, since conversion of the exponent will
affect the fraction part. Given the number $f \cdot 2↑e$ to
be converted to decimal, we may express $2↑e$ in the form $F
\cdot 10↑E$ (usually by means of auxiliar tables), and then
convert $Ff$ to decimal. Alternatively, we can multiply $e$
by log$↓{10} 2 and round this to the nearest integer E$; then
divide $f \cdot 2↑e$ by 10$↑E$ and convert the result. Conversely,
given the number $F \cdot 10↑E$ to be converted to binary, we
may convert $F$ and then multiply it by the floating-point number
10$↑E$ (again commonly obtained by using tables). Obvious techniques
may be used to reduce the maximum size of the auxiliary tables
by using several multiplications and/or divisions, although
this can cause rounding errors to propagate.

\subsectionbegin{E. Multiple-precision conversion}
When converting multiple-precision numbers, it is most convenient
to start by converting blocks of digits, which can be handled
by singe-precision techniques, then to combine them using simple
multiple-precision techniques. For example, suppose that 10$↑n$
is the highest power of 10 less than the computer word size.
Then
$$\quad a) To convert a multiple-precision {\sl integer} from
binary to decimal, divide repeatedly by 10$↑n$ (thus converting
from binary to radix 10$↑n$ by Method (1a)). Single-precision
operations will give the $n$ decimal digits for each place of
the radix-10$↑n$ representation.

b) To convert a multiple-precision {\sl fraction} from binary
to decimal, proceed similarly, multiplying by 10$↑n$ (i.e.,
using Method (2a) with $B = 10↑n)$.

c) To convert a multiple-precision integer from
decimal to binary, convert blocks of $n$ digits first; then
convert from radix 10$↑n$ to binary by using Method (1b).

d) Finally, a multiple-precision fraction may be converted from
decimal to binary by a technique similar to (c), using Method
(2b).

\subsectionbegin{F. History and Bibliography} Radix-conversion
techniques implicitly originated in ancient problems
dealing with weights, measures, and currencies, when a mixed-radix
system was generally involved. Auxiliary tables were usually
used as an aid to conversion. During the seventeenth century
when sexagesimal fractions were being supplanted by decimal
fractions, it was necessary to convert between the two systems
in order to use existing books of tables; a systematic method
to transform fractions from radix 60 to radix 10 and vice versa
was given in the 1667 edition of William Oughtred's {\sl Clavis
Mathematic \xskip }, Chapter 6, Section 18. (This material was
not present in the original 1631 edition of Oughtred's book.)
Conversion rules had been given already by al-Kash|spose 7↓??
of Persia in his {\sl key to Arithmetic} (c. 1414), where Methods
(1a), (1b), and (2a) are clearly explained {\sl [Istoriko-Mat.\
Issled.\ \bf 7} (1954), 126--135], but his work was unknown in
Europe. The 18th century American mathematician Hugh Jones used
the words ``octavation'' and ``decimation'' to describe octal/decimal
conversions, but his methods were not as clever as his terminology.
A. M. Legendre [$Th|spose 1eorie des nombres$ (Paris, 1798),
229] noted that positive integers may be conveniently converted
to binary form by repeatedly dividing by 64.

In 1946, H. H. Goldstine and J. von Neumann gave prominent consideration
to radix conversion in their classic memoir, ``Planning and
coding problems for an electronic computing instrument,'' because
it was necessary to justify the use of binary arithmetic; see
John von Neumann, {\sl Collected Works \bf 5} (New York: Macmillan,
1963), 127--142. Another early discussion of radix conversion
on binary computers was published by F. Koons and S. Lubkin,
{\sl Math.\ Comp.\ \bf 3} (1949), 427--431, where a rather unusual
method is suggested. The first discussion of floating-point
conversion was given somewhat later by F. L. Bauer and K. Samelson
[{\sl Zeit.\ f\"ur angewandte Math.\ und Physik \bf 4} (1953),
312--316].

The following articles may be useful for further reference:
A note by G. T. Lake [{\sl CACM \bf 5} (1962), 468--469] mentions
some hardware techniques for conversion and gives clear examples.
A. H. Stroud and D. Secrest [{\sl Comp.\ J. \bf 6} (1963), 62--66]
have discussed conversion of multiple-precision floating-point
numbers. The conversion of {\sl unnormalized} floating-point
numbers, preserving the amount of ``significance'' implied by
the representation, has been discussed by H. Kanner [{\sl JACM
\bf 12} (1965), 242--246] and by N. Metropolis and R. L. Ashenhurst
[{\sl Math.\ Comp.\ \bf 19} (1965), 435--441]. See also K. Sikdar,
{\sl Sankhy\=a} series B, {\bf 30} (1968), 315--334, and the references
he cites.

|qleft \24skip {\sl }{\:r EXERCISES

\exno 1. [25] Generalize
Method (1b) so that it works with mixed-radix notations, converting
$$$a↓mb↓{m-1} \ldotsm b↓1b↓0 + \cdots + a↓1b↓0 + a↓0\quad$ into\quad
$A↓MB↓{M-1} \ldotsm B↓1B↓0 + \cdots + A↓1B↓0 + A↓0$,
|qctr \9skip where 0 ≤ $a↓j < b↓j, 0 ≤ A↓J < B↓J$ for $0 ≤ j
< m, 0 ≤ J < M$.
|qleft \qquad Give an example of your generalization by manually
converting the quantity ``3 days, 9 hours, 12 minutes, and 37
seconds'' into long tons, hundredweights, stones, pounds, and
ounces. (Let one second equal one ounce. The British system
of weights has 1 stone = 14 pounds, 1 hundredweight = 8 stone,
1 long ton = 20 hundredweight.) In other words, let $b↓0 = 60,
b↓1 = 60, b↓2 = 24, m = 3, B↓0 = 16, B↓1 = 14, B↓2 = 8, B↓3
= 20, M = 4$; the problem is to find $A↓4, \ldotss , A↓0$ in
the proper ranges such that $3b↓2b↓1b↓0 + 2b↓1b↓0 + 12b↓0 +
37 = A↓4B↓3B↓2B↓1B↓0 + A↓3B↓2B↓1B↓0 + A↓2B↓1B↓0 + A↓1B↓0 + A↓0$,
using a systematic method that generalizes Method (1b). (All
arithmetic is to be done in a mixed-radix system.)

\exno 2. [25] Generalize Method (1a) so that it works with mixed-radix
notations, as in exercise 1, and give an example of your generalization
by manually solving the same conversion problem stated in exercise
1.

\exno 3. [25] (D. Taranto.)\xskip The text observes that when fractions
are bing converted, there is in general no obvious way to decide
how many digbits to give in the answer. Design a simple generalization
of Method (2a) that, given two positive radix-$b$ fractions
$u$ and $\cdot$ between 0 and 1, converts $u$ to a radix-$B$
equivalent $U$ that has just enouth places of accuracy to ensure
that $|U - u| < \cdot $. (If $u < \cdot $, we may take $U =
0$, with zero ``places of accuracy.'')

\exno 4. [M21] (a) Prove that every real number having a
terminating {\sl binary} representation also has a terminating
{\sl decimal} representation.\xskip (b) Find a simple condition on
the positive integers $b$ and $B$ that is satisfied if and
only if every real number which has a terminating radix-$b$
representation also has a terminating radix-$B$ representation.

\exno 5. [M20] Show that program (4) would work also if the
instruction ``LDX=10$↑N='' were replaced by ``LDX =$c$='' for
certain other constants $c$.

\exno 6. [30] Discuss using Methods (1a), (1b), (2a), and (2b)
when $b$ or $B$ is -2.

\exno 7. [M18] Given that 0 < $α ≤ x ≤ α + 1/w$ and 0 ≤ $u ≤
w$, prove that \lfloor $ux\rfloor = \lfloor αu\rfloor$ or $\lfloor
αu\rfloor + 1$. Furthermore $\lfloor ux\rfloor = \lfloor αx\rfloor$
exactly, if $u < αw$ and $α↑{-1}$ is an integer.

\exno 8. [24] Write a \MIX\ program analogous to (1) that uses
(5) and includes no division instructions.

\exno 9. [M27] Let $u$ be an integer, 0 ≤ $u < 2↑{34}$. Assume
that the following sequence of operations (equivalent to addition
and binary ``shift-right'' instructions) is performed:
$$
|qctr ⊗v ← \lfloor ${1\over 2}$u\rfloor ,⊗v ← v + \lfloor ${1\over
2}$\rfloor ,⊗v ← v + \lfloor 2$↑{-4}v\rfloor ,\cr
\4skip ⊗v ← v + \lfloor 2↑{-8}v\rfloor ,⊗v ← v + \lfloor 2↑{-16}v\rfloor
,⊗v ← \lfloor {1\over 8}v\rfloor .\cr
\6skip$ Prove that $v = \lfloor u/10\rfloor$ or $\lfloor u/10\rfloor
- 1$.

\exno 10. [22] The text shows how a binary-coded decimal number
can be doubled by using various shifting, extracting, and addition
operations on a binary computer. Give an analogous method that
computes {\sl half} of a binary-coded decimal number (throwing
away the remainder if the number is odd).
|qleft β???????????\quad ??\quad ??\quad ??\quad ??\quad ??\quad
??\quad ??\quad ??\quad ??\quad ??\quad
%folio 411 galley 4a (C) Addison-Wesley 1978	¬
|newtype 58320---Computer Programming\quad (Knuth/Addison-Wesley)\quad
f. 411\quad ch. 4\quad g. 4c
$${\bf 11. \xskip} [{\sl 16}] Convert (57721)$↓8 to decimal$.

\exno 12. [22] Invent a rapid pencil-and-paper method for converting
integers from ternary notation to decimal, and illustrate your
method by converting (1212011210210)$↓3 into decimal. How would
you go fromdecimal to ternary?$

\exno 13. [25] Assume that locations U + 1, U + 2, \ldotss ,
U + $m$ contain a multiple-precision fraction $(.u↓{-1}u↓{-2}
\ldotsm j↓{-m})↓b$, where $b$ is the word size of \MIX: Write
a \MIX\ routine that converts this fraction to decimal notation,
truncating it to 180 decimal digits. The answer should be printed
on two lines, with the digits grouped into 20 blocks of nine
each separated by blanks. (Use the CHAR instruction.)

\exno 14. [M27] (A. Sch\"onhage.)\xskip The text's method of
converting multiple-precision integers requires an execution
time of order $n↑2$ to convert an $n$-digit integer, when $n$
is large. Show that it is possible to convert $n$-digit decimal
integers into binary notation in $O\biglp M(n)$log $n\bigrp$
cycles, where $M(n)$ is an upper bound on the number of cycles
needed to multiply $n$-bit binary numbers that
represents the
best $p$-digit base $b$ floating-point approximation to $u$,
in the sense of Section 4.2.2. Under the assumption that log$↓B
b$ is irrational and that the range of exponents is unlimited,
prove that

|qleft  \9skip $u =$ round$↓b\biglp$ round$↓B(u, P), p\bigrp
$,
|qctr \9skip for all $p$-digit base $b$ floating-point numbers
$u$, if and only if $B↑{p-1} ≥ b↑p$. (In othert words, an ``ideal''
input conversion of $u$ into an independent base $B$, followed
by an ``ideal'' output conversion of this result, will always
yield $u$ again if and only if the intermediate precision $P$
is suitably large, as specified by the formula above.)

\exno 19. [M23] Let the decimal number $u = (u↓7 \ldotsm u↓1u↓0)↓{10}$
be represented as the binary-coded decimal number $U = (u↓7
\ldotsm u↓1u↓0)↓{16}$. Find appropriate constants $c↓i$ and masks
$m↓i$ so that the operation $U ← U - c↓i(U ∧ m↓i)$, repeated
for $i = 1, 2, 3$, will convert $U$ to the binary representation
of $u$, where ``∧'' denotes extraction (i.e., ``logical AND''
on individual bits).

%folio 412 galley 4b (C) Addison-Wesley 1978	-
|qleft \25skip |newtype {\:r 4.5. RATIONAL ARITHMETIC

|qleft }\12skip {\:r 4.5.1. Fractions

|qleft }\6skip In many numerical algorithms it is important
to know that the answer to a problem is exactly ${1\over 3}$,
not a floating-point number that gets printed as ``0.333333574''.
If arithmetic is done on fractions instead of on approximations
to fractions, many computations can be done entirely {\sl without
any accumulated rounding errors.} This results in a comfortable
feeling of security that is often lacking when floating-point
calculations have been made, and it means that the accuracy
of the calculation cannot be improved upon.

When fractional arithmetic is desired, th e numbers can be represented
as pairs of integers, $(u/u↑\prime )$, where $u$ and $u↑\prime$
are relatively prime to each other and $u↑\prime > 0$. The number
zero is represented as (o/1). In this form, $(u/u↑\prime ) =
(v/v↑\prime )$ if and only if $u = v\hfill$ and $u↑\prime =
v↑\prime $.

Multiplication of fractions is, of course, rather
simple; to form $(u/u↑\prime ) \times (v/v↑\prime ) = (w/w↑\prime
)$, we can simply compute $uv$ and $u↑\prime v↑\prime $. The
two products $uv$ and $u↑\prime v↑\prime$ might not be relatively
prime, but if $d =$ gcd$(uv, u↑\prime v↑\prime )$, the desired
answer is $w = uv/d, w↑\prime =u↑\prime v↑\prime /d$. (See exercise
2.) Efficient algorithms to compute the greatest common divisor
are discussed in Section 4.5.2.

Another way to perform the multiplication is to find $d↓1 =$
gcd$(u, v↑\prime )$ and $d↓2 =$ gcd$(u↑\prime , v)$; then the
answer is $w = (u/d↓1)(v/d↓2), w↑\prime = (u↑\prime /d↓2)(v↑\prime
/d↓1)$. (See exercise 3.) This method requires two gcd calculations,
but it is not really slower than the former method; the gcd
process involves a number of iterations essentially proportional
to the logarithm of its inputs, so the total number of iterations
needed to evaluate both $d↓1$ and $d↓2$ is essentially the same
as the number of iterations during the single calculation of
$d$. Furthermore, each iteration in the elevation of $d↓1$ and
$d↓2$ is potentially faster, because comparatively small numbers
are being examined. If $u, u↑\prime , v, v↑\prime$ are single-precision
quantities, this method hs the advantagbe that no double-precision
numbers appear in the calculation unless it is impossible to
represent both of the answers $w$ and $w↑\prime$ in single-precision
form.

Division may be done in a similar manner; see exercise 4.

Addition and subtraction are slightly more complicated. We can
set $(u/u↑\prime ) \pm (v/v↑\prime ) = \biglp (uv↑\prime \pm
u↑\prime v)/u↑\prime v↑\prime \bigrp$ and reduce this fraction
to lowest terms by calculating $d =$ gcd$(uv↑\prime \pm u↑\prime
v, u↑\prime v↑\prime )$ as in the first multiplication method.
But again it is possible to avoid working with such large numbers,
if we start by calculating $d↓1 =$ gcd$(u↑\prime , v↑\prime
)$. If $d↓1 = 1$ then $w = uv↑\prime \pm u↑\prime v$ and $w↑\prime
= ↑\prime v↑\prime$ are the desired answers. (According to THhe???
fraction to lowest terms by calculating $d =$ gcd$(uv↑\prime
\pm u↑\prime v, u↑\prime v↑\prime )$ as in the first multiplication
method. But again it is possible to avoid working with such
large numbers, if we start by calculating $d↓1 =$ gcd$(u↑\prime
, v↑\prime )$. If $d↓1 = 1$ then $w = uv↑\prime \pm u↑\prime
v$ and $w↑\prime = ↑\prime v↑\prime$ are the desired answers.
(According to Theorem 4.5.2D\null, $d↓1$ will be 1 about 61 percent
of the time, if the denominators $u↑\prime$ and $v↑\prime$ are
randomly distributed, so it is wise to single out this case
separately.) If $d↓1 > 1$, then let $t = u(v↑\prime /d↓1)\pm
v(u↑\prime /d↓1)$ and calclate $d↓2 =$ gcd$(t, d↓1)$; finally
the answer is $w = t/d↓2, w↑\prime = (u↑\prime /d↓1)(v↑\prime
/d↓2)$. (Exercise 6 proves that these values of $w$ and $w↑\prime$
are relatively prime to each other.) If single-precision numbers
are being used, this method requires only single-precision operations,
except that $t$ may be a double-precision number or slightly
larger (see exercise 7); since gcd$(t, d↓1) =$ gcd($t$ mod $d↓1,
d↓1)$, the calculation of $d↓2$ does not require double precision.

For example, to compute (7/66) + (17/12), we form $d↓1 =$ gcd(66,
12) = 6; then $t = 7 \cdot 2 + 17 \cdot 11 = 201$, and $d↓2
=$ gcd(201, 6) = 3, so the answer is
$$${201\over 3}$/${66\over 6}$ ${12\over 3}$ = 67/44.

|qleft \9skip \quad Experience with fractional calculations
shows that in many cases the numbers grow to be quite large;
for example, see the remarks following Eq.\ 3.3.2--22. So if
$u$ and $u↑\prime$ are intended to be single-precision numbers
for each fraction $(u/u↑\prime )$, it is important to include
tests for overflow in each of the addition, subtraction, multiplication,
and division subroutines. For numerical problems in which perfect
accuracy is important, a set of subroutines for fractional arithmetic
with {\sl arbitrary} precision allowed in numerator and denominator
is very useful.

The methods of this section also extend to other number fields
besides the rational numbers; for example, we could do arithmetic
on quantities of the form $(u + u↑\prime |newtype |newtype \sqrt{5)u↑{\prime\prime}
$, where $u, u↑\prime , u↑{\prime\prime}$ are integers, gcd$(u, u↑\prime
, u↑{\prime\prime} ) = 1$, and $u↑{\prime\prime} > 0$; or on quantities of the
form $(u + u↑\prime |spose ↑3|newtype |newtype \sqrt{2} + u↑{\prime\prime}
|spose ↑3|newtype |newtype \sqrt{4})/u $, etc.

To help check out subroutines for rational arithmetic, inversion
of matrices with known inverses (e.g., Cauchy matrices, exercises
111111??????\quad To help check out subroutines for rational
arithmetic, inversion of matrices with known inverses (e.g.,
Cauchy matrices, exercise 1.2.3--41) is suggested.

Exacty representation of fractions within a computer was first
discussed in the literature by P. Henrici, {\sl JACM \bf 3}
(1956), 6--9.

|qleft \24skip {\:r EXERCISES

\exno 1. [15] Suggest
a reasonable computational method for comparing two fractions,
to test whether or not $(u/u↑\prime ) < (v/v↑\prime )$.

\exno 2. [M15] Prove that if $d
=$ gcd$(u, v)$ then $u/d$ and $v/d$ are relatively prime.

\exno 3. [M20] Prove that if $u$ and $u↑\prime$ are relatively
prime, and if $v$ and $v↑\prime$ are relatively prime, then
gcd($uv, u↑\prime v↑\prime ) =$ gcd$(u, v↑\prime )$gcd$(u↑\prime
, v)$.

\exno 4. 11] Design a division algorithm for fractions, analogous
to the second multiplication method of the text. (Note that
the sign of $\omega$ must be considered.)

\exno 5. [10] Compute (17/120) + (-27/70) by the method recommended
in the text.

\exno 6. [M23] Show that if $u, u↑\prime$ are relatively prime
and if $v, v↑\prime$ are relatively prime, then gcd($uv↑\prime
+ vu↑\prime , u↑\prime v↑\prime ) = d↓1d↓2$ where $d↓1 =$ gcd$(u↑\prime
, v↑\prime )$ and $d↓2 =$ gcd($d↓1, u(v↑\prime /d↓1) + v(u↑\prime
/d↓1))$. (Hence if $d↓1 = 1$, thenn $εεεεε??????↑\prime /d↓1))$.
(Hence if $d↓1 = 1$, then $uv↑\prime + vu↑\prime$ is relatively
prime to $u↑\prime v↑\prime .)$

\exno 7. [M22] How large can the
absolute value of the quantity $t$ become, in the addition-subtraction
method recommended in the text, if the input numerators and
denominators are less than $N$ in absolute value?

\exno 8. [22] Discuss using $(1/0)$ and $(-1/0)$ as representations
for $∞$ and $-∞$, and/or as representations of overflow.

\exno 9. [M23] If 1 ≤ $u↑\prime
, v↑\prime < 2↑n$, show that \lfloor 2$↑{2n}u/u↑\prime \rfloor
= \lfloor 2↑{2n}v/v↑\prime \rfloor$ implies $u/u↑\prime = v/v↑\prime
$.

\exno 10. [41] Extend the subroutines suggested in exercise
4.3.1--34 so that they deal with ``arbitrary'' rational numbers.

\exno 11. [M23] Consider fractions of the form $(u + u↑\prime
|newtype |newtype \sqrt{5})/u↑{\prime\prime} $, where $u, u↑\prime ,
u↑{\prime\prime}$ are integers, gcd($u, u↑\prime , u↑{\prime\prime} ) = 1$,
and $u↑{\prime\prime} > 0$. Explain how to divide two such fractions
and
%folio 415 galley 1 (C) Addison-Wesley 1978	-
|newtype 58320---Computer Programming\quad (A.-W./Knuth)\quad
f. 415\quad ch. 4\quad g. 1

|qleft \12skip {\:r 4.5.2 The Greatest Common Divisor

|qleft }\6skip If $u$ and $v$ are integers, not both zero, we
say that their {\sl greatest common divisor}, gcd($u, v)$, is
the largest integer that evenly divides both $u$ and $v$. This
definition makes sense, because if $u ≠ 0$ then no integer greater
than |$u|$ can evenly divide $u$, but the integer 1 does divide
both $u$ and $v$; hence there must be a largest integer that
divides them both. When $u$ and $v$ are both zero, every integer
evenly divides zero, so the above definition does not apply;
it is convenient to set
|qleft
 |qctr \6skip \hfill gcd(0, 0) ⊗= 0.\eqno (1)\cr
\6skip \quad The definitions just given obviously imply that
$$\hfill gcd($u, v) ⊗=$ gcd($v, u),\eqno (2)\cr
\6skip \hfill$ gcd($u, v) ⊗=$ gcd(-u, v),\eqno (3)\cr
\6skip \hfill$ gcd($u, 0) ⊗=|u|.\eqno (4)\cr
\6skip$ \quad In the previous section, we reduced the problem
of expressing a rational number in ``lowest terms'' to the problem
of finding the greatest common divisor of its numerator and
denominator. Other applications of the greatest common divisor
have been mentioned for example in Sections 3.2.1.2, 3.3.3,
4.3.2, 4.3.3. So the concept of the greatest common divisor
is important and worthy of serious study.
|qleft ??????\quad The {\sl least common multiple} of two integers
$u$ and $v$, written lcm($u, v)$, is a related idea that is
also important. It is defined to be the smallest positive integer
that is a multiple of (i.e., evenly divisible by) both $u$
and $v$; and lem(0, 0) = 0. The classical method for teaching
children how to add fractions $u/u↑\prime + v/v↑\prime$ is to
train them to find the ``least common denominator,'' which is
lcm($u↑\prime , v↑\prime )$.

According to the ``fundamental theorem of arithmetic''
(proved in exercise 1.2.4--21), each positive integer $u$ can
be expressed in the form
$$$u = 2↑u|infsup 23↑u|infsup 35↑u|infsup 57↑u|infsup 711↑u|infsup
1|infsup 1 \ldots = \prod ↓{p$ prime} $p↑u|infsup p,\eqno (5)$
|qctr \6skip where the exponents $u↓2, u↓3, . . $. are uniquely
determined nonnegative integers, and where all but a finite
number of the exponents are zero. From this canonical factorization
of a positive integer, it is easy to discover one way to compute
the greatest common divisor of $u$ and $v:$ By (2), (3), and
(4) we may assume that $u$ and $v$ are positive integers, and
if both of them have been canonically factored into primes,
we have
$$gcd($u, v) |tab = \prod ↓{p$ prime} $p↑{$min(u$↓p,v↓p)},\eqno
(6)⊗\6skip \hfill lcm(u, v) ⊗= \prod ↓{p$ prime} $p↑{$max($u↓p,v↓p)}.\eqno
(7)\cr$
Thus, for example, the greatest common divisor of $u = 7000
= 2↑3 \cdot 5↑3 \cdot 7$ and $v = 4400 = 2↑4 \cdot 5↑2 \cdot
11$ is 2↑{min(3,4)}5↑{min(3,2)}7↑{min(1,0)}11↑{min(0,1)} = 2$↑3
\cdot 5↑2 = 200. The least common multiple of the same two numbers
is 2↑4 \cdot 5↑3 \cdot 7 \cdot 11 = 154000$.

From formulas (6) and (7) we can easily prove a
number of basic identities concerning the gcd and the lcm:
$$
|qleft \hfill gcd($u, v)w ⊗=$gcd($uw, vw),⊗$if⊗$w ≥ 0;\eqno
(8)\cr
\6skip \hfill$ lcm($u, v)w ⊗=$ lcm($uw, vw),⊗$if⊗$w ≥ 0;\eqno
(9)\cr
\6skip \hfill u \cdot v ⊗=$ gcd($u, v) \cdot$ lcm($u, v),⊗$if⊗$u,
v ≥ 0;\eqno (10)\cr
\6skip \hfill$ gcd\biglp lcm($u, v)$, lcm($u, w)\bigrp ⊗=$ lcm\biglp
u, gcd($v, w)\bigrp ;\eqno (11)\cr
\6skip \6skip$ The latter two formulas are ``distinctive laws''
analogous to the familiar identity $uv + uw = u(v + w)$. Equation
(10) reduces the calculation of gcd($u, v)$ to the calculation
of lcm($u, v)$, and conversely.

\subsectionbegin{Euclid's algorithm} Although Eq.\
(6) is useful for theoretical purposes, it is generally no help
for calculating a greatest common divisor in practice, because
it requires that we first determine the factorization of $u$
and $v$. There is no known method for finding the prime factors
of an integer very rapidly (see Section 4.5.4). But fortunately
there is an efficient way to calculate the greates common divisor
of two integers without factoring them, and, in fact, such a
method was discovered over 2250 years ago; this is ``Euclid's
algorithm,'' which we have already examined in Sections 1.1
and 1.2.1.

Euclid's algorithm is found in Book 7, Propositions
1 and 2 of his {\sl Elements (c. 300 |newtype} B|newtype .|newtype
C|newtype .), but it probably wasn't his own invention. Scholars
believe that the method was known up to 200 years earlier, at
least in its subtractive form, and it was almost certainly known
to Eudoxus $c. 375 |newtype$ B|newtype .|newtype C|newtype .
[Cf.\ K. von Fritz, {\sl Ann.\ Math.\ (2) \bf 46} (1945), 242--264.]
We might call it the granddaddy of all algorithms, because it
is the oldest nontrivial algorithm that has survived to the
present day. (The chief rival for this honor is perhaps the
ancient Egyptian method for multiplication, which was based
on doubling and adding, and which forms the basis for efficient
calculation of $n$th powers as explained in Section 4.6.3. But
the Egyptian manuscripts merely give examples that are not
completely systematic, and that were certainly not stated systematically;
the Egyptian method is therefore not quite deserving of the
name ``algorithm.'' Several ancient Babylonian methods, for
doing such things as solving certain sets of quadratic equations
in two variables, are also known. Genuine algorithms are involved
in this case, not just special solutions to the equations for
certain input parameters; even though the Babylonians invariably
presented each method in conjunction with an example worked
with special parameters, they regularly explained the general
procedure in the accompanying text. [See D. E. Knuth, {\sl CACM \bf 15}
(1972), 671677.] Many of these Babylonian algorithms predate
Euclid by 1500 years, and they are the earliest known instances
of written procedures for mathematics. They do not have the
stature of Euclid's algorithm, however, since they do not involve
iteration and since they have been superseded by modern algebraic
methods.)

Since Euclid's algorithm is there fore important
for historical as well as practical reasons, let us |newtype
58320---Co
%folio 417 galley 2 (C) Addison-Wesley 1978	-
mputer Programming\quad (A.-W./Knuth)\quad f. 417\quad ch. 4\quad
g. 2

|qleft \12skip |newtype |indent {\bf Proposition. \xskip} {\sl
Given two positive integers, find their greatest common divisor.}

|qleft \12skip Let $A, C$ be the two given positive integers;
it is required to find their greatest common divisor. If $C$
divides $A$, then $C$ is a common divisor of $C$ and $A$, since
it also divides itself. And it clearly is in fact the greatest,
since no greater number than $C$ will divide $C$.

|qleft \12skip But if $C$ does not divide $A$, then continually
subtract the lesser of the numbers $A, C$ from the greater,
until some number is left that divides the previous one. This
will eventually happen, for if unity is left, it will divide
the previous number.

|qleft \12skip Now let $E$ be the positive remainder of $A$
divided by $C$; let $F$ be the positive remainder of $C$ divided
by $E$; and let $F$ be a divisor of $E$. Since $F$ divides $E$
and $E$ divides $C - F, F$ also divides $C - F$; but it also
divides itself, so it divides $C$. And $C$ divides $A - E$;
therefore $F$ also divides $A - E$. But it also divides $E$;
therefore it divides $A$. Hence it is a common divisor of $A$
and $C$.

|qleft \12skip I now claim that it is also the greatest. For
if $F$ is not the greatest common divisor of $A$ and $C$, some
larger number will divide them both. Let such a number be $G$.

|qleft \12skip Now since $G$ divides $C$ while $C$ divides $A
- E, G$ divides $A - E. G$ also divides the whole of $A$, so
it divides the remainder $E$. But $E$ divides $C - F$; therefore
$G$ also divides $C - F. G$ also divides the whole of $C$, so
it divides the remainder $F$; that is, a greater number divides
a smaller one. This is impossible.

|qleft \12skip Therefore no number greater than $F$ will divide
$A$ and $C$, and so $F$ is their greatest common divisor.

|qleft \12skip {\bf Corollary. \xskip} This argument makes it
evident that any number dividing two numbers divides their greatest
common divisor. {\sl Q.E.D.}

|qleft \12skip |cancelindent |newtype Note. \xskip Euclid's
statements have been simplified here in one nontrivial respect:
Greek mathematicians did not regard unity as a ``divisor'' of
another positive integer. Two positive integers were either
both equal to uhity, or they were relatively prime, or they
had a greawt?????????hey were relatively prime, or they had
a greatest common divisor. In fact, unity was not even considered
to be a ``number,'' and zero was of course nonexistent. These
rather awkward conventions made it necessary for Euclid to duplicate
much of his discussion, and he gave two separate propositions
that each are essentially like the one appearing here.

In his discussion, Euclid first suggests subtracting
the smaller of the two current numbers from the larger repeatedly,
until we get two numbers in which one is a multiple of another;
but in the proof he really relies on taking the remainder of
one number divided by another (and since he has no concept of
zero, he cannot speak of the remainder when one number divides
the other). So it is reasonable to say that he imagines each
{\sl division} (not the individual subtractions) as a single
step of the algorithm, and hence an ``authentic'' rendition
of his algorithm can be phrased as follows:

\algbegin Algorithm E (Original Euclidean algorithm).
Given two integers {\sl A and C} greater than unity, this algorithm
finds their greatest common divisor.

\algstep E1. [$A$ divisible by $C$?] If $C$ divides
$A$, the algorithm terminates with $C$ as the answer.

\algstep E2. [Replace $A$ by remainder.] If $A$ mod $C$ is equal
to unit, the given numbers were relatively prime, so the algorithm
terminates. Otherwise replace the pair of values $(A, C)$ by
$(C, A$ mod $C)$ and return to step E1.
|qleft |cancelindent \12skip \quad The ``proof'' Euclid gave,
which is quoted above, is especially interesting because it
is, of course, not really a proof at all! He verifies the result
of the algorithm only if step E1 is performed once or thrice.
Surely he must have realized that step E1 could take place more
than three times, although he made no mention of such a possibility.
Not having the notion of a proof by mathematical induction,
he could only give a proof for a finite number of cases. (In
fact, he often proved only the case $n = 3$ of a theorem that
he wanted to establish for general $n.)$ Although Euclid is
justly famous for the great advances he made in the art of logical
deduction, techniques for giving valid proofs by indukc?????????e
hot discovered until many centuries later, and the crucial ideas
for proving the validity of {\sl algorithms} are only now becoming
really clear. (See Section 1.2.1 for a complete proof of Euclid's
algorithm, together with a short discussion of general proof
procedures for algorithms.)

It is worth noting that this algorithm for finding
the greatest common divisor was chosen by Euclid to be the very
first step in his development of the theory of numbers. The
same order of presentation is still in use today in modern textbooks.
Euclid also gave (Proposition 34) a method to find the least
common multiple of two integers $u$ and $v$, namely to divide
$u$ by gcd($u, v)$ and to multiply by $n$; this is equivalent
to Eq.\ (10).

If we avoid Euclid's bias against the numbers 0
and 1, we can reformulate Algorithm E in the following way:

\algbegin Algorithm A (Modern Euclidean algorithm).
Given nonnegative integers $u$ and $v$, this algorithm finds
their greatest common divisor. ({\sl Note{\sl : }}The greatest
common divisor of {\sl arbitrary} integers $u$ and $v$ may be
obtained by applying this algorithm to $|u|$ and $|v|$, because
of Eqs.\ (2) and (3)|newtype )

\algstep A1. [$v = 0$?] If $v = 0$, the
algorithm terminates with $u$ as the answer.

\algstep A2. [Take $u \mod v$.] Set
$r ← u$ mod $v, u ← v, v ← r$, and return to A1. The operations
of this step decrease the value of $v$, but leave gcd($u, v)$
unchanged.)

|qleft \12skip \quad For example, we may calculate gcd(40902,
24140) as follows:
|qleft |cancelindent \12skip
 |qctr \hfill gcd(40902, 24140) ⊗= gcd(24140, 16762) = gcd(16762,
7378)\cr
\4skip ⊗= gcd(7378, 2006) = gcd(2006, 1360) = gcd(1360, 646)\cr
\4skip ⊗= gcd(646, 68) = gcd(68, 34) = gcd(34, 0) = 34.\cr
\6skip \quad A proof that Algorithm A is valid follows readily
from Eq.\ (4) and the fact that
$$gcd($u, v) =$ gcd($v, u - qv),\eqno (13)$
|qctr \6skip if $q$ is any integer. Equation (13) holds because
any common divisor of $u$ and $v$ is a divisor of both $v$ and
$u - qv$, and, conversely, any common divisor of $v$ and $u
- qv$ must divide both $u$ and $v$.

The following \MIX\ program illustrates the fact
that Algorithm A can easily be implemented on a computer:

\algbegin Program A (Euclid's algorithm). Assume that
$u$ and $v$ are single-precision, nonnegative integers, stored
respectively in locations U and V; this program puts gcd($u,
v)$ into rA.

|qleft \6skip |newtype |tab \qquad |tab 333 |tab \qquad \quad
|tab \qquad \qquad \quad |tab \qquad \qquad \qquad \qquad \quad
|tab |zfa ⊗⊗LDX⊗U⊗1⊗rX ← $u.\cr$
⊗JMP⊗2F⊗1⊗\cr
1H⊗STX⊗V⊗$T⊗v ←$ rX.\cr
⊗SRAX⊗5⊗??U0}|newtype 58320---Computer
%folio 419 galley 3 (C) Addison-Wesley 1978	-
\subsectionbegin{A binary method} Since Euclid's
patriarchal algorithm has been used for so many centuries, it
is a rather surprising fact that it may not always be the best
method for finding the greatest common divisor after all. A
quite different gcd algorithm, which is primarily suited to
binary arithmetic, was discovered by J. Stein in 1961 [see {\sl
J. Comp.\ Phys.\ \bf 1} (1967), 397--405.] This new algorithm requires
no division instruction; it i based sslely on the operations of subtraction,
testing whether a number is even or odd, and shifting the binary
representation of an even number to the right (halving).

The binary gcd algorithm is based on four simple facts about positive
integers $u$ and $v$:

\yskip\textindent{a)}If $u$ and $v$ are both even, then $\gcd(u,
v) = 2\gcd(u/2, v/2)$. [See Eq.\ (8).]

\textindent{b)}If $u$ is even and $v$ is odd, then $\gcd(u, v) =
\gcd(u/2, v)$. [See Eq.\ (6).]

\textindent{c)}As in Euclid's algorithm, $\gcd(u, v) = \gcd(u - v,
v)$. [See Eqs.\ (13), (2).]

\textindent{d)}If $u$ and $v$ are both odd, then $u - v$ is even,
and $|u - v| < \max(u, v)$.

\yskip\noindent These facts immediately suggest the following
algorithm:

\algbegin Algorithm B (Binary gcd algorithm). Given
positive integers $u$ and $v$, this algorithm finds their greatest
common divisor.

\algstep B1. [Find power of 2.] Set $k ←
0$, and then repeatedly set $k ← k + 1$, $u ← u/2$, $v ← v/2$ zero
or more times until $u$ and $v$ are not both even.

\algstep B2. [Initialize.] (Now the original values of $u$ and
$v$ have been divided by $2↑k$, and at least one of their present
values is odd.) If $u$ is odd, set $t ← -v$ and go to B4. Otherwise
set $t ← u$.

\algstep B3. [Halve $t$.] (At this point, $t$ is even,
and nonzero.) Set $t ← t/2$.

\algstep B4. [Is $t$ even?] If $t$ is even, go back to B3.

\algstep B5. [Reset $\max(u, v)$.] If $t > 0$, set $u ← t$;
otherwise set $v ← -t$. (The larger of $u$ and $v$ has been
replaced by $|t|$, except perhaps during the first time this
step is performed.)

\algstep B6. [Subtract.] Set $t ← u - v$. If $t ≠ 0$, go back
to B3. Otherwise the algorithm terminates with $u \cdot 2↑k$
as the output.\quad\blackslug

\topinsert{\vskip 48mm
\ctrline{\caption Fig.\ 9. Binary algorithm for the greatest common divisor.}}

\yyskip As an example of Algorithm
B\null, let us consider $u = 40902$, $v = 24140$, the same numbers
we have used to try out Euclid's algorithm. Step B1 sets $k
← 1$, $u ← 20451$, $v ← 12070$. Then $t$ is set to $-12070$, and replaced
by $-6035$; $v$ is replaced by 6035, and the computation proceeds
as follows:
$$\vjust{\halign{\hfill#⊗\qquad\hfill#⊗\qquad#\hfill\cr
$u$\hfill⊗$v$\hfill⊗\hfill$t$\cr
\noalign{\vskip 3pt}
20451⊗6035⊗$ +14416$, $+7208$, $+3604$, $+1802$, $+901$;\cr
901⊗6035⊗$-5134$, $-2567$;\cr
901⊗2567⊗$-1666$, $-833$;\cr
901⊗833⊗$+68$, $+34$, $+17$;\cr
17⊗833⊗$-816$, $-408$, $-204$, $-102$, $-51$;\cr
17⊗51⊗$-34$, $-17$;\cr
17⊗17⊗0.\cr}}$$
The answer is $17 \cdot 2↑1 = 34$. A few more iterations
were necessary here than we needed with Algorithm A\null, but each
iteration was somewhat simpler since no division steps were
used.

A \MIX\ program for Algorithm B requires just a little
more code than for Algorithm A\null. In order to make such a program
fairly typical of a binary computer representation of Algorithm
B\null, let us assume that \MIX\ is extended to include the following
operators:

\yyskip\noindent$\bullet$ \.{SLB} (shift left AX binary).\xskip
C$\null = 6$; F$\null = 6$. \par\penalty1000\noindent
The contents of registers A and X are ``shifted
left'' M binary places; that is, $|\rAX| ← |2↑{\char'115}\rAX|\mod B↑{10}$,
where $B$ is the byte size. (As with all \MIX\ shift commands,
the signs of rA and rX are not affected.)

\yyskip\noindent$\bullet$ \.{SBR} (shift
right AX binary).\xskip
C$\null = 6$; F$\null = 7$.\par\penalty1000\noindent
The contents of registers A
and X are ``shifted right'' M binary places; that is, $|\rAX|
← \lfloor |rAX|/2↑{\char'115}\rfloor $.

\yyskip\noindent$\bullet$ \.{JAE}, \.{JAO} (jump
A even, jump A odd).\xskip C$\null = 40$; F$\null = 6$, 7, respectively.\par
\penalty1000\noindent
A \.{JMP} occurs if rA is even or odd, respectively.

\yyskip\noindent$\bullet$ \.{JXE}, \.{JXO} (jump
X even, jump X odd).\xskip C$\null = 47$; F$\null = 6$, 7, respectively.\par
\penalty1000\noindent
Analogous to \.{JAE}, \.{JAO}.

\algbegin Program B (Binary gcd algorithm).
Assume that $u$ and $v$ are single-precision positive integers,
stored respectively in locations \.U and \.V; this program uses
Algorithm B to put $\gcd(u, v)$ into rA. Register assignments:
$t ≡ \rA$, $k ≡ \rI1$.
%folio 421 galley 4 WARNING: Much tape unreadable! (C) Addison-Wesley 1978	-
{\yyskip\tabskip 25pt \mixfive{\!
01⊗ABS⊗EQU⊗1:5\cr
02⊗B1⊗ENT1⊗0⊗1⊗\understep{B1. Find }{\sl p\hskip-3pt}\understep{\hskip 3pt
ower of 2.}\cr
03⊗⊗LDX⊗U⊗1⊗$\rX ← u$.\cr
04⊗⊗LDAN⊗V⊗1⊗$\rA ← -v$.\cr
05⊗⊗JMP⊗1F⊗1\cr
\\06⊗2H⊗SRB⊗1⊗A⊗Halve rA, rX.\cr
07⊗⊗INC1⊗1⊗A⊗$k ← k + 1$.\cr
08⊗⊗STX⊗U⊗A⊗$u ← u/2$.\cr
09⊗⊗STA⊗V(ABS)⊗A⊗$v ← v/2$.\cr
\\10⊗1H⊗JXO⊗B4⊗1+A⊗To B4 with $t←-v$ if $u$ is odd.\cr
\\11⊗B2⊗JAE⊗2B⊗B+A⊗\understep{B2. Initialize.}\cr
\\12⊗⊗LDA⊗U⊗B⊗$t←u$.\cr
\\13⊗B3⊗SRB⊗1⊗D⊗\understep{B3. Halve $t$.}\cr
\\14⊗B4⊗JAE⊗B3⊗1-B+D⊗\understep{B4. Is $t$ even?}\cr
\\15⊗B5⊗JAN⊗1F⊗C⊗\understep{B5. Reset $\max(u,v)$.}\cr
16⊗⊗STA⊗U⊗E⊗If $t>0$, set $u←t$.\cr
17⊗⊗SUB⊗V⊗E⊗$t←u-v$.\cr
18⊗⊗JMP⊗2F⊗E\cr
\\19⊗1H⊗STA⊗V(ABS)⊗C-E⊗If $t<0$, set $v←-t$.\cr
20⊗B6⊗ADD⊗U⊗C-E⊗\understep{B6. Subtract.}\cr
\\21⊗2H⊗JANZ⊗B3⊗C⊗To B3 if $t≠0$.\cr
\\22⊗⊗LDA⊗U⊗1⊗$\rA←u$.\cr
23⊗⊗ENTX⊗0⊗1⊗$\rX←0$.\cr
24⊗⊗SLB⊗0,1⊗1⊗$\rA←2↑k\cdot\rA$.\quad\blackslug\cr}}

The running time of this program is
$$9A+2B+6C+3D+E+13$$
units, where $A=k$, $B=1$ if $t←u$ in step B2 (otherwise $B=0$), $C$ is the
number of subtraction steps, $D$ is the number of halvings in step B3, and
$E$ is the number of times $t>0$ in step B5. Calculations discussed later in
this section imply that we may take $A={1\over3}$, $B={1\over3}$, $C=0.71n-0.5$,
$D=1.41n-2.7$, $E=0.35n-0.4$ as average values for these quantities, assuming
random inputs $u$ and $v$ in the range $1≤u,v2↑n$. The total running time is
therefore about $8.8n+5$ cycles, compared to about $11.1n+7$ for Program A under
the same assumptions. The worst possible running time for $u$ and $v$ in this
range occurs when $A=0$, $B=0$, $C=n$, $D=2n-2$, $E=1$; this amounts to
$12n+8$ cycles. (The corresponding value for Program A is $26.8n+19$.)

Thus the greater speed of the iterations in Program B, due to the simplicity of
the operations, compensates for the greater number of iterations required.
We have found that the binary algorithm is about 20 percent faster than
Euclid's algorithm on the \MIX\ computer. Of course, the situation may be
different on other computers, and in any event both programs are quite
efficient; but it appears that not even a procedure as venerable as Euclid's
algorithm can withstand progress.

V. C. Harris [{\sl Fibonacci Quarterly \bf8} (1970), 102--103] has suggested an
interesting cross between Euclid's algorithm and the binary algorithm. If
$u$ and $v$ are odd, with $u≥v>0$, we can always write $u=qv\pm r$ where
$0≤r<v$ and $r$ is even; if $r≠0$ we set $r←r/2$ until $r$ is odd, then set
$u←v$, $v←r$ and repeat the process. In subsequent ierations, $q≥3$.

\subsectionbegin{Extensions} We can extend the methods
used to calculate $\gcd(u, v)$ in order to solve some slightly
more difficult problems. For example, assume that we want to
compute the greatest common divisor of $n$ integers $u↓1$, $u↓2$,
$\ldotss$, $u↓n$.

One way to calculate $\gcd(u↓1, u↓2, \ldotss , u↓n)$,
assuming that the $u$'s are all nonnegative, is to extend Euclid's
algorithm in the following way: If all $u↓j$ are zero, the greatest
common divisor is taken to be zero; otherwise if only one $u↓j$
is nonzero, it is the greatest common divisor; otherwise replace
$u↓k$ by $u↓k \mod u↓j$ for all $k ≠ j$, where $u↓j$ is the
minimum of the nonzero $u$'s.

The algorithm sketched in the preceding paragraph
is a natural generalization of Euclid's method, and it can be
justified in a similar manner. But there is a simpler method
available, based on the easily verified identity
$$\gcd(u↓1, u↓2, \ldotss , u↓n) =\gcd\biglp u↓1,\gcd(u↓2,
\ldotss , u↓n)\bigrp .\eqno (14)$$
To calculate $\gcd(u↓1, u↓2, \ldotss , u↓n)$, we
may proceed as follows:

\yskip\hang\textindent{\bf D1.}Set $d ← u↓n$, $j ← n - 1$.

\yskip\hang\textindent{\bf D2.}If $d≠1$ and $j>0$, set $d←\gcd(u↓j, d)$
and $j ← j - 1$ and repeat this step. Otherwise
$d = \gcd(u↓1, \ldotss , u↓n)$.

\yskip\noindent This method reduces the calculation
of $\gcd(u↓1, \ldotss , u↓n)$ to repeated calculations of the
greatest common divisor of two numbers at a time. It makes use
of the fact that $\gcd(u↓1, \ldotss , u↓j, 1) = 1$; and this will
be helpful, since we will already have $\gcd(u↓{n-1}, u↓n) =
1$ over 60 percent of the time if $u↓{n-1}$ and $u↓n$ are chosen
at random. In most cases, the value of $d$ will decrease rapidly during
in the first few stages of the calculation, and this will make the
remainder of the computation quite fast. Here Euclid's algorithm
has an advantage over Algorithm B\null, in that its running time
is primarily governed by the value of $\min(u, v)$, while the
running time for Algorithm B is primarily governed by $\max(u,
v)$; it would be reasonable to perform one iteration of Euclid's
algorithm, replacing $u$ by $u \mod v$ if $u$ is much larger
than $v$, and then to continue with Algorithm B.

The assertion that $\gcd(u↓{n-1}, u↓n)$ will be
equal to unity more than 60 percent of the time for random inputs
is a consequence of the following well-known result of number
theory:

\algbegin Theorem D ({\rm G. Lejeune Dirichlet, {\sl Abh.\ K\"oniglich
Preuss.\ Akad.\ Wiss.}\ (1849), 69--83}). {\sl If $u$ and $v$ are integers
chosen at random, the probability that $\gcd(u, v) = 1$ is $6/π↑2$.}

\yyskip A precise formulation of this theorem, which carefully
defines what is meant here by being ``chosen at random,'' appears
in exercise 10 with a rigorous proof. Let us content ourselves
here with a heuristic argument that shows why the theorem is plausible.
%folio 423 galley 5 (C) Addison-Wesley 1978	-
If we assume, without proof, the existence
of a well-defined probability $p$ that $\gcd(u, v)$ equals unity,
then we can determine the probability that $\gcd(u, v) = d$ for
any positive integer $d$; for the latter event will happen only
when $u$ is a multiple of $d$, $v$ is a multiple of $d$, and $\gcd(u/d,
v/d) = 1$. Thus the probability that $\gcd(u, v) = d$ is equal
to $1/d$ times $1/d$ times $p$, namely $p/d↑2$. Now let us sum
these probabilities over all possible values of $d$; we should get
$$1 = \sum ↓{d≥1} p/d↑2 = p(1 + {1\over 4} + {1\over 9} + {1\over
16} +\cdotss ).$$
Since the sum $1 + {1\over 4} + {1\over 9}+
\cdots = H↑{(2)}↓{\!\!∞}$ is equal to $π↑2/6$ (cf.\ Section 1.2.7),
we must have $p=6/π↑2$ in order to make this
equation come out right.\quad\blackslug

\yyskip Euclid's algorithm can be extended
in another important way: We can calculate integers $u↑\prime$
and $v↑\prime$ such that
$$uu↑\prime + vv↑\prime = \gcd(u, v)\eqno (15)$$
at the same time $\gcd(u, v)$ is being calculated.
This extension of Euclid's algorithm can be described conveniently
in vector notation:

\algbegin Algorithm X (Extended Euclid's algorithm).
Given nonnegative integers $u$ and $v$, this algorithm determines
a vector $(u↓1, u↓2, u↓3)$ such that $uu↓1 + vu↓2 = u↓3 = \gcd(u,
v)$. The computation makes use of auxiliary vectors $(v↓1, v↓2,
v↓3)$, $(t↓1, t↓2, t↓3)$; all vectors are manipulated in such
a way that the relations
$$ut↓1 + vt↓2 = t↓3,\qquad uu↓1 + vu↓2 = u↓3,\qquad uv↓1 +
vv↓2 = v↓3\eqno (16)$$
hold throughout the calculation.

\algstep X1. [Initialize.] Set $(u↓1, u↓2, u↓3) ←
(1, 0, u)$, $(v↓1, v↓2, v↓3) ← (0, 1, v)$.

\algstep X2. [Is $v↓3 = 0$?] If $v↓3 = 0$, the algorithm terminates.

\algstep X3. [Divide, subtract.] Set $q ← \lfloor u↓3/v↓3\rfloor
$, and then set
$$\baselineskip15pt
\cpile{(t↓1, t↓2, t↓3) ← (u↓1, u↓2, u↓3) - (v↓1, v↓2, v↓3)q,\cr
(u↓1, u↓2, u↓3) ← (v↓1, v↓2, v↓3),\qquad (v↓1,
v↓2, v↓3) ← (t↓1, t↓2, t↓3).\cr}$$
Return to step X2.\quad\blackslug

\yyskip For example, let $u = 40902$, $v = 24140$.
At step X2 we have
$$\vjust{\tabskip20pt\halign{\hfill#⊗\hfill$#$⊗\hfill$#$⊗\hfill#⊗\hfill$#$⊗\hfill
$#$⊗\hfill#\cr
$q$⊗u↓1⊗u↓2⊗$u↓3$\hfill⊗v↓1⊗v↓2⊗$v↓3$\hfill\cr
---⊗1⊗0⊗40902⊗0⊗1⊗24140\cr
1⊗0⊗1⊗24140⊗1⊗-1⊗16762\cr
1⊗1⊗-1⊗16762⊗-1⊗2⊗7378\cr
2⊗-1⊗2⊗7378⊗3⊗-5⊗2006\cr
3⊗3⊗-5⊗2006⊗-10⊗17⊗1360\cr
1⊗-10⊗17⊗1360⊗13⊗-22⊗646\cr
2⊗13⊗-22⊗646⊗-36⊗61⊗68\cr
9⊗-36⊗61⊗68⊗337⊗-571⊗34\cr
2⊗337⊗-571⊗34⊗-710⊗1203⊗0\cr}}$$
Therefore the solution is $337 \cdot 40902 - 571 \cdot
24140 = 34 = \gcd(40902, 24140)$.

\yakip The validity of Algorithm X follows from
(16) and the fact that the algorithm is identical to Algorithm
A with respect to its manipulation of $u↓3$ and $v↓3$. A detailed
proof of Algorithm X is discussed in Section 1.2.1. Gordon H.
Bradley has observed that we can avoid a good deal of the calculation
in Algorithm X by suppressing $u↓1$, $v↓1$, and $t↓1$; then $u↓1$
can be determined afterwards using the relation $uu↓1 + vu↓2
= u↓3$. (It is interesting to note that this modified algorithm would
be harder to verify; we often find that the simplest way to prove an
``optimized'' algorithm correct is to show that it is equivalent to an
unoptimized algorithm whose correctness is easier to establish.)

Exercise 14 shows that the values of $|u↓1|$, $|u↓2|$,
$|v↓1|$, $|v↓2|$ remain bounded by the size of the inputs $u$ and
$v$.

Algorithm B\null, which computes the greatest common
divisor properties of binary notation, can be extended in a similar
(but more complicated) way; see exercise 35. For some
instructive extensions to Algorithm X\null, see exercises 18 and
19 in Section 4.6.1.

\yskip The ideas underlying Euclid's algorithm
can also be applied to find a {\sl general solution in integers}
of any set of linear equations with integer coefficients. For
example, suppose that we want to find all integers $w$, $x$, $y$,
$z$ that satisfy the two equations
$$\vjust{\tabskip0pt plus 100pt\vskip-9pt\baselineskip16pt
\halign to size{$\hfill#$\tabskip0pt⊗$\null#⊗$\null\hfill$⊗$\null#$⊗$\null#\hfill
\tabskip0pt plus 100pt⊗\hfill#\cr
10w⊗+ 3x⊗+ 3y⊗+ 8z⊗= 1,⊗(17)\cr
 6w⊗- 7x⊗    ⊗- 5z⊗= 2.⊗(18)\cr}}$$
We can introduce a new variable
$$\lfloor 10/3\rfloor w + \lfloor 3/3\rfloor x + \lfloor 3/3\rfloor
y + \lfloor 8/3\rfloor z = 3w + x + y + 2z = t↓1,$$
and use it to eliminate $y$; Eq.\ (17) becomes
$$(10 \mod 3)w + (3 \mod 3)x + 3t↓1 + (8 \mod 3)z = w + 3t↓1
+ 2z = 1,\eqno (19)$$
and Eq.\ (18) remains unchanged. The new equation
(19) may be used to eliminate $w$, and (18) becomes
$$6(1 - 3t↓1 - 2z) - 7x - 5z = 2;$$
that is,
$$7x + 18t↓1 + 17z = 4.\eqno (20)$$
Now as before we introduce a new variable
$$x + 2t↓1 + 2z = t↓2$$
and eliminate $x$ from (20):
$$7t↓2 + 4t↓1 + 3z = 4.\eqno (21)$$
Another new variable can be introduced in the same
fashion, in order to eliminate the variable $z$, which has the
smallest coefficient:
$$2t↓2 + t↓1 + z = t↓3.$$
Eliminating $z$ from (21) yields
$$t↓2 + t↓1 + 3t↓3 = 4,\eqno (22)$$
and this equation, finally, can be used to eliminate
$t↓2$. We are left with two independent variables, $t↓1$ and
$t↓3$; substituting back for the original variables, we obtain
the general solution
$$\baselineskip14pt
\eqalign{w ⊗= \quad 17 -\9 5t↓1 - 14t↓3,\cr
x ⊗=\quad 20 - \9 5t↓1- 17t↓3,\cr
y ⊗= -55 + 19t↓1 + 45t↓3,\cr
z ⊗=\9 -8 +\quad t↓1 +\9 7t↓3.\cr}\eqno(23)$$
In other words, all integer solutions
$(w, x, y, z)$ to the original equations (17), (18) are obtained
from (23) by letting $t↓1$ and $t↓3$ independently run through
all integers.
%folio 426 galley 6 (C) Addison-Wesley 1978	-

The general method that has just been illustrated
is based on the following procedure: Find a nonzero coefficient
$c$ of smallest absolute value in the system of equations. Suppose
that this coefficient appears in an equation having the form
$$cx↓0+c↓x↓1+\cdots+c↓kx↓k=d;\eqno(24)$$
we may assume for simplicity that $c > 0$. If
$c = 1$, use this equation to eliminate the variable $x↓0$ from
the other equations remaining in the system; then repeat the
procedure on the remaining equations. (If no more equations
remain, the computation stops, and a general solution in terms
of the variables not yet eliminated has essentially been obtained.)
If $c > 1$, then if $c↓1 \mod c =\cdots = c↓k \mod
c = 0$ we must have $d \mod c = 0$, otherwise there is no
integer solution; divide both sides of (24) by $c$ and eliminate
$x↓0$ as in the case $c = 1$. Finally, if $c > 1$ and not all
of $c↓1\mod c$, $\ldotss$, $c↓k\mod c$ are zero, then introduce
a new variable
$$\lfloor c/c\rfloor x↓0 + \lfloor c↓1/c\rfloor
x↓1 +\cdots + \lfloor c↓k/c\rfloor x↓k = t;\eqno(25)$$
eliminate the variable $x↓0$ from the other equations,
in favor of $t$, and replace the original equation (24) by
$$ct + (c↓1 \mod c)x↓1 +\cdots+ (c↓k\mod c)x↓k= d.\eqno (26)$$
$\biglp$Cf.\ (19) and (21) in the above example.$\bigrp$

This process must terminate, since each step either
reduces the number of equations or the size of the smallest
nonzero coefficient in the system. A study of the above procedure
will reveal its intimate connection with Euclid's algorithm.
The method is a comparatively simple means of solving linear
equations when the variables are required to take on integer
values only. It isn't the best available method for this problem,
however; substantial refinements are possible, but beyond the scope of
this book.

\subsectionbegin{High-precision calculation} If $u$
and $v$ are very large integers, requiring a multiple-precision
representation, the binary method (Algorithm B) is a simple
and reasonably efficient means of calculating their greatest
common divisor.

By contrast, Euclid's algorithm seems much less
attractive, since step A2 requires a multiple-precision division
of $u$ by $v$. But this difficulty is not really as bad as it
seems, since we will prove in Section 4.5.3 that the quotient
$\lfloor u/v\rfloor$ is almost always very small; for example,
assuming random inputs, the quotient $\lfloor u/v\rfloor$ will
be less than 1000 approximately 99.856 percent of the time.
Therefore it is almost always possible to find $\lfloor u/v\rfloor$
and $(u\mod v)$ using single-precision calculations, together
with the comparatively simple operation of calculating $u -
qv$ where $q$ is a single-precision number.

A significant improvement in the speed of Euclid's
algorithm when high-precision numbers are involved can be achieved
by using a method due to D. H. Lehmer [{\sl AMM \bf 45}
(1938), 227--233]. Working only with the leading digits of large
numbers, it is possible to do most of the calculations with
single-precision arithmetic, and to make a substantial reduction
in the number of multiple-precision operations involved.

Lehmer's method can be illustrated on the eight-digit
numbers $u = 27182818$, $v = 10000000$, assuming that we are using
a machine with only four-digit words. Let $u↑\prime = 2718$,
$v↑\prime = 1001$, $u↑{\prime\prime} = 2719$, $v↑{\prime\prime} = 1000$;
then $u↑\prime/v↑\prime$ and $u↑{\prime\prime} /v↑{\prime\prime}$
are approximations to $u/v$, with
$$u↑\prime /v↑\prime < u/v < u↑{\prime\prime} /v↑{\prime\prime} .\eqno (27)$$
The ratio $u/v$ determines the sequence of quotients
obtained in Euclid's algorithm; if we carry out Euclid's algorithm
simultaneously on the single-precision values $(u↑\prime , v↑\prime
)$ and $(u↑{\prime\prime} , v↑{\prime\prime} )$ until we get a different quotient,
it is not difficult to see that the same sequence of quotients
would have appeared to this point if we had worked with the
multiple-precision numbers $(u, v)$. Thus, consider what happens
when Euclid's algorithm is applied to $(u↑\prime , v↑\prime
)$ and to $(u↑{\prime\prime} , v↑{\prime\prime} )$:
$$\vjust{\halign{\hfill#⊗\qquad\hfill#⊗\qquad\hfill#⊗\hskip80pt\hfill#⊗\qquad
\hfill#⊗\qquad\hfill#\cr
$u↑\prime$\hfill⊗$v↑\prime$\hfill⊗$q↑\prime$\hfill⊗$u↑{\prime\prime}$\hfill
⊗$v↑{\prime\prime}$\hfill⊗$q↑{\prime\prime}$\hfill\cr
\noalign{\vskip 3pt}
2718⊗1001⊗2⊗2719⊗1000⊗2\cr
1001⊗716⊗1⊗10000⊗719⊗1\cr
716⊗285⊗2⊗719⊗281⊗2\cr
285⊗146⊗1⊗281⊗157⊗1\cr
146⊗139⊗1⊗157⊗124⊗3\cr
139⊗7⊗19⊗124⊗33⊗3\cr}}$$
After six steps we find that $q↑\prime ≠ q↑{\prime\prime} $,
so the single-precision calculations are suspended; we have
gained the knowledge that the calculation would have proceeded
as follows if we had been working with the original multiple-precision
numbers:
$$\vjust{\halign{$\ctr{#}$⊗\qquad$\ctr{#}$⊗\qquad$\ctr{#}$\cr
u⊗v⊗q\cr
\noalign{\vskip 3pt}
u↓0⊗v↓0⊗2\cr
v↓0⊗u↓0- 2v↓0⊗1\cr
u↓0- 2v↓0⊗ -u↓0 + 3v↓0⊗2\cr
-u↓0 + 3v↓0⊗ 3u↓0 - 8v↓0⊗1\cr
3u↓0 -8v↓0⊗ -4u↓0 +11v↓0⊗1\cr
-4u↓0 + 11v↓0⊗ 7u↓0 ⊗- 19v↓0⊗?\cr}}\eqno(28)$$
(The next quotient lies somewhere between 3 and 19.)
No matter how many digits are in $u$ and $v$, the first five
steps of Euclid's algorithm would be the same as (28), so long
as (27) holds. We can therefore avoid the multiple-precision
operations of the first five steps, and replace them all by
a multiple-precision calculation of $-4u↓0 + 11v↓0$ and $7u↓0
- 19v↓0$. In this case we obtain $u = 1268728$, $v = 279726$;
the calculation can now proceed with $u↑\prime = 1268$, $v↑\prime
= 280$, $u↑{\prime\prime} = 1269$, $v↑{\prime\prime} = 279$, etc. With a larger
accumulator, more steps could be done by single-precision calculations;
our example showed that only five cycles of Euclid's algorithm
were combined into one multiple step, but with (say) a word
size of 10 digits we could do about twelve cycles at a time.
(Results proved in Section 4.5.3 imply that the number of multiple-precision
cycles that can be replaced at each iteration is essentially
proportional to the logarithm of the word size used in the single-precision
calculations.)

Lehmer's method can be formulated as follows:

\algbegin Algorithm L (Euclid's algorithm for large numbers).
Let $u, v$ be nonnegative integers, with $u ≥ v$, represented
in multiple precision. This algorithm computes the greatest
common divisor of $u$ and $v$, making use of auxiliary single-precision
$p$-digit variables $u$, $v$, $A$, $B$, $C$, $D$, $T$, and $q$, and auxiliary
multiple-precision variables $t$ and $w$.

\algstep L1. [Initialize.] If $v$ is small
enough to be represented as a single-precision value, calculate
$\gcd(u, v)$ by Algorithm A and terminate the computation. Otherwise,
let $\A u$ be the $p$ leading digits of $u$, and let $\A v$
be the $p$ corresponding digits of $v$; in other words,
if radix-$b$ notation is being used, $\A u←\lfloor u/b↑k\rfloor$,
\A v← \lfloor v/b↑k\rfloor $, where $k$ is as small
as possible consistent with the condition $\A u < b↑p$.

\hangindent 19pt after 0 Set $A ← 1$, $B ← 0$, $C ← 0$, $D ← 1$. (These
variables represent the coefficients in (28), where
$$u = Au↓0 + Bv↓0,\qquad v = Cu↓0 + Dv↓0,\eqno (29)$$
in the equivalent actions of Algorithm
A on multiprecision numbers. We also have
$$u↑\prime = \A u + B,\qquad v↑\prime =
\A v + D,\qquad u↑{\prime\prime} = \A u + A,\qquad v↑{\prime\prime} = 
\A v + C\eqno (30)$$
in terms of the notation in the example worked above.)

\algstep L2. [Test quotient.] Set $q ← \lfloor (\A u +
A)/(\A v + C)\rfloor $. If $q ≠ \lfloor (\A u + B)/(
\A v + D)\rfloor $, go to step L4. (This step tests if $q↑\prime
≠ q↑{\prime\prime}$ in the notation of the above example. Note that
single-precision overflow can occur in special circumstances
during the computation in this step, but only when $\A u
= b↑p - 1$ and $A = 1$ or when $\A v = b↑p - 1$ and $D
= 1$; the conditions
$$\baselineskip15pt
\eqalign{0⊗≤\A u + A ≤ b↑p,\cr 0⊗≤\A u + B < b↑p,\cr}\qquad
\eqalign{0⊗≤\A v + C < b↑p,\cr 0⊗≤\A v + D ≤ b↑p\cr}\eqno(31)$$
will always hold, because of (30). It
is possible to have $\A v + C = 0$ or $\A v + D =
0$, but not both simultaneously; therefore division by zero
in this step is taken to mean ``Go directly to L4.'')

%folio 429 galley 7 (C) Addison-Wesley 1978	-
\algstep L3. [Emulate Euclid.] Set $T ← A - qC$, $A ← C$, $C ← T$,
$T ← B - qD$, $B ← D$, $D ← T$, $T ← \A u - q\A v$,
$\A u ← \A v$, $\A v ← T$, and go back to step L2. $\biglp$These
single-precision calculations are the equivalent of multiple-precision
operations, as in (28), under the conventions of (29).$\bigrp$

\algstep L4. [Multiprecision step.] If $B = 0$, set $t ← u
\mod v$, $u ← v$, $v ← t$, using multiple-precision division. (This
happens only if the single-precision operations cannot simulate
any of the multiprecision ones. It implies that Euclid's algorithm
requires a very large quotient, and this is an extremely rare
occurrence.) Otherwise, set $t ← Au$, $t ← t + Bv$, $w ← Cu$, $w ←
w + Dv$, $u ← t$, $v ← w$ (using straightforward multiprecision
operations). Go back to step L1.\quad\blackslug

\yyskip The values of $A$, $B$, $C$, $D$
remain as single-precision numbers throughout this calculation,
because of (31).

Algorithm L requires a somewhat more complicated program
than Algorithm B\null, but with large numbers it will be faster on
many computers. Algorithm B can, however, be speeded up in a
similar way (see exercise 34), to the point where it continues
to win. Algorithm L has the advantage that it can be extended,
as in Algorithm X (see exercise 17); furthermore, it determines
the sequence of quotients obtained in Euclid's algorithm, and
this yields the regular continued fraction expansion of a real
number (see exercise 4.5.3--18).

\subsectionbegin{Analysis of the binary algorithm}
Let us conclude this section by studying the running time of
Algorithm B\null, in order to justify the formulas stated earlier.

An exact determination of the behavior of Algorithm
B appears to be exceedingly difficult to derive, but we can
begin to study it by means of an approximate model of its behavior.
Suppose that $u$ and $v$ are odd numbers, with $u > v$ and
$$\lfloor\lg u\rfloor = m,\qquad \lfloor\lg v\rfloor
= n.\eqno(32)$$
$\biglp$Thus, $u$ is an $(m + 1)$-bit number, and
$v$ is an $(n + 1)$-bit number.$\bigrp$ Algorithm B forms $u
- v$ and shifts this quantity right until obtaining an odd number
$u↑\prime$ that replaces $u$. Under random conditions, we would
expect to have
$u↑\prime = (u - v)/2$
about one-half of the time,
$u↑\prime = (u - v)/4$
about one-fourth of the time,
$u↑\prime = (u - v)/8$
about one-eighth of the time, and so on. We have
$$\lfloor\lg u↑\prime \rfloor = m - k - r,\eqno (33)$$
where $k$ is the number of places that $u - v$
is shifted right, and where $r$ is $\lfloor\lg u\rfloor -
\lfloor\lg (u - v)\rfloor$, the number of bits lost at
the left during the subtraction of $v$ from $u$. Note that $r
≤ 1$ when $m ≥ n + 2$, and $r ≥ 1$ when $m = n$. For simplicity,
we will assume that $r = 0$ when $m ≠ n$ and that $r = 1$ when
$m = n$, although this lower bound tends to make $u↑\prime$
seem larger than it usually is.

The approximate model we shall use to study Algorithm
B is based solely on the values $m = \lfloor\lg u\rfloor$
and $n = \lfloor\lg v\rfloor$ throughout the course of the
algorithm, not on the actual values of $u$ and $v$. Let us call
this approximation a {\sl lattice-point model}, since we will
say that we are ``at the point $(m, n)$'' when $\lfloor\lg
u\rfloor = m$ and $\lfloor\lg v\rfloor = n$. From point $(m,
n)$ the algorithm takes us to $(m↑\prime , n)$ if $u > v$, or
to $(m, n↑\prime )$ if $u < v$, or terminates if $u = v$. For
example, the calculation starting with $u = 20451$, $v = 6035$
that is tabulated after Algorithm B begins at the point $(14,
12)$, then goes to $(9, 12)$, $(9, 11)$, $(9, 9)$, $ 4, 9)$, $(4, 5)$,
$(4, 4)$, and terminates. In line with the comments of the preceding
paragraph, we will make the following assumptions about the
probability that we reach a given point just after point $(m,
n)$:

\yyskip
\hjust to size{\vjust{\hjust to 14.5pc{\hfill Case 1, $m>n$.\hfill}
\vskip 3pt
\tabskip 0pt plus 100pt
\halign to 14.5pc{$\ctr{#}$⊗$\ctr{#}$\cr
\hjust{Next point}⊗\hjust{Probability}\cr
\noalign{\vskip 2pt}
(m-1,n)⊗1\over2\cr
(m-2,n)⊗1\over4\cr
\cdots⊗\cdots\cr
(1,n)⊗{1\over2}↑{m-1}\cr
(0,n)⊗{1\over2}↑{m-1}\cr}}\!
\vjust{\hjust to 14.5pc{\hfill Case 2, $m<n$.\hfill}
\vskip 3pt
\tabskip 0pt plus 100pt
\halign to 14.5pc{$\ctr{#}$⊗$\ctr{#}$\cr
\hjust{Next point}⊗\hjust{Probability}\cr
\noalign{\vskip 2pt}
(m,n-1)⊗1\over2\cr
(m,n-2)⊗1\over4\cr
\cdots⊗\cdots\cr
(m,1)⊗{1\over2}↑{n-1}\cr
(m,0)⊗{1\over2}↑{n-1}\cr}}}
\hjust to size{\hfill Case 3, $m=n>0$.\hfill}
\vskip 3pt
\tabskip 0pt plus 100pt
\halign to size{$\ctr{#}$⊗$\ctr{#}$\cr
\hjust{Next point}⊗\hjust{Probability}\cr
\noalign{\vskip 2pt}
(m-2,n),\,(m,n-2)⊗{1\over4},\,{1\over4}\cr
(m-3,n),\,(m,n-3)⊗{1\over8},\,{1\over8}\cr
\cdots⊗\cdots\cr
(0,n),\,(m,0)⊗{1\over2}↑m,\,{1\over2}↑m\cr
\hjust{terminate}⊗{1\over2}↑{m-1}\cr}}
For example, from points $(5, 3)$ the lattice-point model
would go to points $(4, 3)$, $(3, 3)$, $(2, 3)$, $(1, 3)$, $(0, 3)$ with
the respective probabilities ${1\over 2}$, ${1\over 4}$, ${1\over
8}$, ${1\over 16}$, ${1\over 16}$; from $(4, 4)$ it would go to
$(2, 4)$, $(1, 4)$, $(0, 4)$, $(4, 2)$, $(4, 1)$, $(4, 0)$, or would terminate,
with the respective probabilities ${1\over 4}$, ${1\over 8}$,
${1\over 16}$, ${1\over 4}$, $1\over8$, ${1\over 16}$, ${1\over
8}$. When $m = n = 0$, the formulas above do not apply; the
algorithm always terminates in such a case, since $m = n = 0$
implies that $u = v = 1$.

The detailed calculations of exercise 18 show that
this lattice-point model is somewhat pessimistic. In fact, when
$m > 3$ the actual probability that Algorithm B goes from $(m,
m)$ to one of the two points $(m - 2, m)$ or $(m, m - 2)$ is
equal to ${1\over 8}$, although we have assumed the value ${1\over
2}$; the algorithm actually goes from $(m, m)$ to $(m - 3, m)$
or $(m, m - 3)$ with probability ${7\over 32}$, not ${1\over
4}$; it actually goes from $(m + 1, m)$ to $(m, m)$ with probability
${1\over 8}$, not ${1\over 2}$. The probabilities in the model
are nearly correct when $|m - n|$ is large, but when $|m - n|$
is small the model predicts slower convergence than is actually
obtained. In spite of the fact that our model is not a completely
faithful representation of Algorithm B\null, it has one great virtue,
namely that it can be completely analyzed! Furthermore, empirical
experiments with Algorithm B show that the behavior predicted
by the lattice-point model is analogous to the true behavior.

An analysis of the lattice-point model
can be carried out by solving the following rather complicated
set of double recurrence relations:
$$\vcenter{\baselineskip22pt\halign{$\rt{#}$⊗$\dispstyle\null#,$\hfill\qquad
⊗if $#\hfill$\cr
A↓{mm}⊗= a + {1\over 2}A↓{m(m-2)} +\cdots
+ {1\over 2↑{m-1}} A↓{m0} + {b\over 2↑{m-1}}⊗m≥ 1;\cr
A↓{mn}⊗= c + {1\over 2}A↓{(m-1)n} +\cdots
+ {1\over 2↑{m-1}} A↓{1n} + {1\over 2↑{m-1}} A↓{0n}⊗m > n ≥ 0;\cr
A↓{mn} = A↓{nm}⊗n > m ≥ 0.\cr}}\eqno(34)$$
The problem is to solve for $A↓{mn}$ in terms
of $m$, $n$, and the parameters $a$, $b$, $c$ and $A↓{00}$. This is
an interesting set of recurrence equations, which have an interesting
solution.

First we observe that if $0 ≤ n < m$,
$$\eqalign{A↓{(m+1)n} ⊗= c + \sum ↓{1≤k≤m}2↑{-k}A↓{(m+1-k)n}
+ 2↑{-m}A↓{0n}\cr
⊗= c + {1\over 2}A↓{mn} + \sum ↓{1≤k<m}2↑{-k-1}A↓{(m-k)n}
+ 2↑{-m}A↓{0n}\cr
⊗= c + {1\over 2}A↓{mn} + {1\over 2}(A↓{mn} - c)\cr
\noalign{\vskip 4pt}
⊗= {c\over 2} + A↓{mn}.\cr}$$
Hence $A↓{(m+k)n}={1\over2}ck+A↓{mn}$, by induction on $k$. In particular, 
since $A↓{10}=c+A↓{00}, we have
$$\textstyle A↓{m0}={1\over2}c(m+1)+A↓{00},\qquad m>0.\eqno(35)$$
%folio 432 galley 8 (C) Addison-Wesley 1978	-
Now let $A↓m = A↓{mm}$. If $m > 0$, we have
$$\eqalign{A↓{(m+1)m}⊗ = c + \sum ↓{1≤k≤m+1} 2↑{-k}A↓{(m+1-k)m}
+ 2↑{-m-1}A↓{0m}\cr
 ⊗= c + {1\over 2}A↓{mm} + \sum ↓{1≤k≤m}
\biglp 2↑{-k-1}(A↓{(m-k)(m+1)} - c/2)\bigrp + 2↑{-m-1}A↓{0m}\cr
 ⊗\textstyle= c + {1\over 2}A↓m + {1\over 2}(A↓{(m+1)(m+1)} -
a - 2↑{-m}b) - {1\over 4}c(1 - 2↑{-m})\cr
\noalign{\vskip 2pt}
⊗\textstyle\hskip 65mm\null+ 2↑{-m-1}\biglp{1\over2}c(m + 1) + A↓{00}\bigrp\cr
\noalign{\vskip 5pt}
⊗\textstyle= {1\over 2}(A↓m + A↓{m+1}) + {3\over 4}c
- {1\over 2}a + 2↑{-m-1}(c - b + A↓{00}) + m2↑{-m-2}c.\cr}\eqno(36)$$
Similar maneuvering, as shown in exercise 19,
establishes the relation
$$\textstyle A↓{n+1} = {3\over 4}A↓n + {1\over 4}A↓{n-1} + α + 2↑{-n-1}β
+ (n + 2)2↑{-n-1}\gamma ,\quad n ≥ 2,\eqno (37)$$
where
$$\textstyle α = {1\over 4}a + {7\over 8}c,\qquad β = A↓{00} - b - {3\over
2}c,\qquad\hjust{and}\qquad \gamma = {1\over 2}c.$$

Thus the double recurrence (34) can be transformed
into the single recurrence relation in (37). Use of the generating
function $G(z) = A↓0 + A↓1z + A↓2z↑2 +\cdots$ now
transforms (37) into the equation
$${\textstyle(1 - {3\over 4}z - {1\over 4}z↑2)}G(z) = a↓0 + a↓1z +
a↓2z↑2 + {α\over 1 - z} + {β\over 1 - z/2} + {\gamma \over (1
- z/2)↑2},\eqno(38)$$
where $a↓0$, $a↓1$, and $a↓2$ are constants that
can be determined by the values of $A↓0$, $A↓1$, and $A↓2$. Since
$1 - {3\over 4}z + {1\over 4}z↑2 = (1 + {1\over 4}z)(1 - z)$,
we can express $G(z)$ by the method of partial fractions in
the form
$$G(z) = b↓0 + b↓1z + {b↓2\over (1 - z)↑2} + {b↓3\over 1 -
z} + {b↓4\over (1 - z/2)↑2} + {b↓5\over 1 - z/2} + {b↓6\over
1 + z/4} .$$
Tedious manipulations produce the values of these
constants $b↓0$, $\ldotss$, $b↓6$, and therefore the coefficients
of $G(z)$ are determined. We finally obtain the solution
$$\baselineskip15pt
\eqalignno{A↓{nn} ⊗\textstyle= n({1\over 5}a + {7\over 10}c) + ({16\over
25}a + {2\over 5}b - {23\over 50}c + {3\over 5}A↓{00})\cr
⊗\quad\null + 2↑{-n}(-{1\over 3}cn + {2\over 3}b - {1\over 9}c
- {2\over 3}A↓{00})\cr
⊗\quad\null + (-{1\over 4})↑n(-{16\over 25}a - {16\over 15}b +
{16\over 225}c + {16\over 15}A↓{00}) + {1\over 2}c\delta ↓{n0};\cr
\noalign{\vskip 3pt}
A↓{mn} ⊗\textstyle= {1\over2}mc + n({1\over 5}a + {1\over 5}c)
+ ({6\over 25}a + {2\over 5}b + {7\over 50}c + {3\over 5}A↓{00})
+ 2↑{-n}({1\over 3}c)\cr
⊗\quad\null + (-{1\over 4})↑n(-{6\over 25}a - {2\over 5}b + {2\over
75}c + {2\over 5}A↓{00}),\qquad m > n.⊗(39)\cr}$$

With these elaborate calculations done, we can
readily determine the behavior of the lattice-point model. Assume
that the inputs $u$ and $v$ to the algorithm are odd, and let
$m = \lfloor\lg↓2 u\rfloor$, $n = \lfloor\lg v\rfloor $.
The average number of subtraction cycles, namely, the quantity
$C$ in the analysis of Program B\null, the number of times step B6
is executed, is obtained by setting $a = 1$, $b = 0$, $c = 1$, $A↓{00}
= 1$ in the recurrence (34). By (39) we see that (for $m ≥ n$)
the lattice model predicts
$$\textstyle C = {1\over 2}m + {2\over 5}n + {49\over 50} - {1\over
5}\delta ↓{mn}\eqno (40)$$
subtraction cycles, plus terms that rapidly go
to zero as $n$ approaches infinity.

The average number of times that $\gcd(u, v) = 1$ is
obtained by setting $a = b = c = 0$, $A↓{00} = 1$; this gives
the probability that $u$ and $v$ are relatively prime, approximately
${3\over 5}$. Actually, since $u$ and $v$ are assumed to be
odd, they should be relatively prime with probability $8/π↑2$
(see exercise 13), so this reflects the degree of inaccuracy
of the lattice-point model.

The average number of times that a path from $(m, n)$
goes through one of the ``diagonal'' points $(m↑\prime , m↑\prime
)$ for some $m↑\prime ≥ 1$ is obtained by setting $a = 1$, $b
= c = A↓{00} = 0$ in (34); so we find that this quantity is
approximately
$${1\over 5}n + {6\over 25} + {2\over 5}\delta ↓{mn},\qquad
\hjust{when }m ≥ n.$$
Now we can determine the average number of shifts,
the number of times step B3 is performed. (This is the quantity
$D$ in Program B\null.) In any execution of Algorithm B\null, with $u$
and $v$ both odd, the corresponding path in the lattice model
must satisfy the relation
$$\hjust{number of shifts} + \hjust{number of diagonal points} + 2\lfloor\lg
\gcd(u, v)\rfloor = m + n,$$
since we are assuming that $r$ in (33) is always
either 0 or 1. The average value of $\lfloor\lg gcd(u, v)\rfloor$
predicted by the lattice-point model is approximately ${4\over
5}$ (see exercise 20). Therefore we have, for the total number
of shifts,
$$\eqalign{D⊗\textstyle= m + n - ({1\over 5}n + {6\over 25} + {2\over 5}\delta
↓{mn}) -{4\over 5}\cr
\noalign{\vskip2pt}
⊗=\textstyle m + {4\over 5}n - {46\over 25}- {2\over 5}\delta ↓{mn},\cr}\eqno (41)$$
when $m ≥ n$, plus terms that decrease rapidly to zero
for large $n$.

To summarize the most important facts we have derived
from the lattice-point model, we have shown that if $u$ and
$v$ are odd and if $\lfloor\lg u\rfloor = m$, $\lfloor\lg
v\rfloor = n$, then the quantities $C$ and $D$ that are the
critical factors in the running time of Program B will have
average values given by
$$\textstyle C = {1\over 2}m + {2\over 5}n + O(1),\qquad D = m + {4\over
5}n + O(1),\qquad m ≥ n.\eqno (42)$$
But the model that we have used to derive (42)
is only a pessimistic approximation to the true behavior; Table
1 compares the true average values of $C$, computed by actually
running Algorithm B with all possible inputs, to the values
predicted by the lattice-point model, for small $m$ and $n$.
The lattice model is completely accurate when $m$ or $n$ is
zero, but it tends to be 
less accurate when $|m - n|$ is small and $\min(m, n)$ is
large. When $m = n = 9$, the lattice-point model gives $C =
8.78$, compared to the true value $C = 7.58$.

\topinsert{\tablehead{Table 1}
\yskip
\ctrline{NUMBER OF SUBTRACTIONS IN ALGORITHM B}
\yskip
\baselineskip 0pt
\halign to size{#\hfill\quad⊗#\tabskip 100pt⊗#\tabskip0pt⊗#⊗\quad\hfill#\cr
\lower 3pt\null⊗\hfill$n$\hfill⊗⊗\hill$n$\hfill\cr
⊗\leaders\hrule\hfill⊗⊗\leaders\hrule\hfill\cr
\lower 3pt\vjust to 15pt{}\!
⊗\hfill0\hfill\hfill1\hfill\hfill2\hfill\hfill3\hfill\hfill4\hfill\hfill5\hfill
⊗⊗\hfill0\hfill\hfill1\hfill\hfill2\hfill\hfill3\hfill\hfill4\hfill\hfill5\hfill\cr
\\0⊗1.00\quad2.00\quad2.50\quad3.00\quad3.50\quad4.00⊗\!
⊗1.00\quad2.00\quad2.50\quad3.00\quad3.50\quad4.00⊗0\cr
\\1⊗2.00\quad1.00\quad2.50\quad3.00\quad3.50\quad4.00⊗\!
⊗2.00\quad1.00\quad3.00\quad2.75\quad3.63\quad3.94⊗1\cr
\\2⊗2.50\quad2.50\quad2.25\quad3.38\quad3.88\quad4.38⊗\!
⊗2.50\quad3.00\quad2.00\quad3.50\quad3.88\quad4.25⊗2\cr
\\3⊗3.00\quad3.00\quad3.38\quad3.25\quad4.22\quad4.72⊗\!
⊗3.00\quad2.75\quad3.50\quad2.88\quad4.13\quad4.34⊗3\cr
\\4⊗3.50\quad3.50\quad3.88\quad4.22\quad4.25\quad5.10⊗\!
⊗3.50\quad3.63\quad3.88\quad4.13\quad3.94\quad4.80⊗4\cr
\\5⊗4.00\quad4.00\quad4.38\quad4.72\quad5.10\quad5.19⊗\!
⊗4.00\quad3.94\quad4.25\quad4.34\quad4.80\quad4.60⊗5\cr
$m$\lower 7pt\vjust to 15pt{}⊗\hfill Predicted by model\hfill
⊗⊗\hfill Actual average values\hfill⊗$m$\cr
⊗\leaders\hrule\hfill⊗⊗\leaders\hrule\hfill\cr
by model⊗⊗Actual average values⊗$m$\cr}}
%folio 434 galley 9 (C) Addison-Wesley 1978	-
Empirical tests of Algorithm B with several
million random inputs and with various values of $m, n$ in the
range $29 ≤ m, n ≤ 37$ indicate that the actual average behavior
of the algorithm is given by
$$\baselineskip 15pt
\rpile{C\approx\null\cr D\approx\null\cr}
\rpile{\textstyle{1\over 2}m\cr m\cr}
\lpile{\null + 0.203n\cr \null+0.41n\cr}
\lpile{\null + 1.9 - 0.4(0.6)↑{m-n},\cr \null - 0.5 - 0.7(0.6)↑{m-n},\cr}
\qquad m ≥ n.\eqno (43)$$
These experiments showed a rather small standard
deviation from the observed average values. The coefficients
${1\over 2}$ and 1 of $m$ in (42) and (43) can be verified rigorously
without using the lattice-point approximation (see exercise
21); so the error in the lattice-point model is apparently in
the coefficient of $n$, which is too high.

The above calculations have been made under the
assumption that $u$ and $v$ are odd and in the ranges $2↑m <
u < 2↑{m+1}$, $2↑n < v < 2↑{n+1}$. If we say instead that $u$
and $v$ are to be {\sl any} integers, independently and uniformly
distributed over the ranges
$$1 ≤ u < 2↑N,\qquad 1 ≤ v < 2↑N,$$
then we can calculate the average values of $C$
and $D$ from the data already given; in fact, if $C↓{mn}$ denotes
the average value of $C$ under our earlier assumptions, exercise
22 shows that we have
$$\twoline{\hskip-10pt(2↑N - 1)↑2C = N↑2C↓{00}+2N\sum↓{1≤n≤N}(N-n)2↑{n-1}
C↓{n0}}{0pt}{\null+2\hskip-2pt\sum↓{1≤n<m≤N}\hskip-2pt(N-m)(N-n)2↑{m+n-2}C↓{mn}
+\hskip-1pt\sum↓{1≤n≤N}\hskip-2pt(N-n)↑22↑{2n-2}C↓{nn}.\quad(44)\hskip-10pt}$$
The same formula holds for $D$ in terms of $D↓{mn}$.
If the indicated sums are carried out using the approximations
in (43), we obtain
$$C \approx 0.70N + 0(1),\qquad D \approx 1.41N + O(1).$$
(See exercise 23.) This agrees perfectly with the
results of further empirical tests, made on several million
random inputs for $N ≤ 30$; the latter tests show that we may
take
$$C = 0.70N - 0.5,\qquad D = 1.41N - 2.7\eqno (45)$$
as good estimates of the values, given this distribution
of the inputs $u$ and $v$.

Richard Brent has discovered a continuous model
that accounts for the leading terms in (45). Let us assume
that $u$ and $v$ are large, and that the number of shifts per
iteration has the value $d$ with probability exactly $2↑d$.
If we let $X = u/v$, the effect of steps B3--B5 is to replace
$X$ by $(X - 1)/2↑d$ if $X > 1$, or by $2↑d/(X↑{-1} - 1)$ if $X <
1$. The random variable $X$ has a limiting distribution that
makes it possible to analyze the average value of the ratio
by which $\max(u, v)$ decreases at each iteration; see exercise
25. Numerical calculations show that this maximum decreases
by $β = 0.705971246102$ bits per iteration; the agreement with
experiment is so good that Brent's constant $β$ must be the
true value of the number ``0.70'' in (45), and we should replace
0.203 by 0.206 in (43). [See {\sl Algorithms and Complexity}, ed.\ by
J. F. Traub (New York: Academic Press, 1976), 321--355.]

This completes our analysis of the average values
of $C$ and $D$. The other three quantities appearing in the
running time of Algorithm B are rather easily analyzed; see
exercises 6, 7, and 8.

Thus we know approximately how Algorithm B behaves
on the average. Let us now consider a ``worst case'' analysis:
What values of $u$ and $v$ are in some sense the hardest to handle?
If we assume as before that
$$\lfloor\lg u\rfloor = m\qquad\hjust{and}\qquad \lfloor
\lg v\rfloor = n,$$
let us try to find $(u, v)$ that make the algorithm
run most slowly. In view of the fact that the subtractions take
somewhat longer than the shifts, when the auxiliary bookkeeping
is considered, this question may be phrased by asking for the inputs
$u$ and $v$ that require most subtractions. The answer is somewhat
surprising; the maximum value of $C$ is exactly
$$\max(m, n) + 1,\eqno (46)$$
although the lattice-point model would predict
that substantially higher values of $C$ are possible (see exercise
26). The derivation of the worst case (46) is quite interesting,
so it has been left as an amusing problem for the reader to
work out by himself (see exercises 27, 28).

\exbegin{EXERCISES}

\exno 1. [M21] How can
(8), (9), (10), (11), and (12) be derived easily from (6) and
(7)?

\exno 2. [M22] Given that $u$ divides $v↓1v↓2 \ldotsm v↓n$, prove
that $u$ divides
$$gcd($u, v↓1)$ gcd($u, v↓2) \ldotsm\gcd($u, v↓n)$.

\exno 3. [M23] Show that the number
of ordered pairs of positive integers $(u, v)$ such that lcm($u,
v) = n$ is the number of divisors of $n↑2$.

\exno 4. [M21] Given positive integers
$u$ and $v$, show that there are divisors $u↑\prime$ of $u$
and $v↑\prime$ of $v$ such that gcd($u↑\prime , v↑\prime ) =
1$ and $u↑\prime v↑\prime =$ lcm($u, v)$.

\exno 5. [M26] Invent an algorithm
(analogous to Algorithm B) for calculating the greatest common
divisor of two integers based on their {\sl balanced ternary}
representation. Demonstrate your algorithm by applying it to
the calculation of gcd(40902,24140).

\exno 6. [M22] Given that $u$ and $v$ are random positive integers,
find the mean and standard deviation of the quantity $A$ (the
number of right shifts on both $u$ and $v$ during the preparatory
phase), which enters into the timing of Program B.

\exno 7. [M20] Analyze the quantity $B$ that enters into the
timing of Program B.

\exno 8. [M24] Show that in Program B\null, the average value of
$E$ is approximately equal to ${1\over 2}$$C$$↓{ave}, where
C$$↓{ave} is the average value of C$.

\exno 9. [18] Using Algorithm B
and hand calculation, find gcd(31408,2718). Also find integers
$m$ and $n$ such that 31408$m + 2718n =$ gcd(31408,2718), using
Algorithm X.

\exno 10. [HM24] Let $q↓n$ be the number of ordered pairs of
integers $(u, v)$ such that $1 ≤ u, v ≤ n$ and gcd($u, v) =
1$. The object of this exercise is to prove that lim$↓{n→∞}
q↓n/n↑2 = 6/π↑2$, thereby establishing Theorem D\null.
|qleft \qquad a) Use the principle of inclusion and exclusion
(Section 1.3.3) to show that
$$$q↓n = n↑2 - \sum ↓{p↓1} Ln/p↓1\rfloor ↑2 + \sum ↓{p↓1<p↓2}
\lfloor n/p↓1p↓2\rfloor ↑2 - \cdotss$,
|qctr \6skip where the sums are taken over all {\sl prime} numbers
$p↓i$.
|qleft \qquad b) The $M|spose 4obius function \mu (n)$ is defined
by the rules $\mu (1) = 1, \mu (p↓1p↓2 \ldotsm p↓r) = (-1)↑r$
if $p↓1, p↓2, \ldotss , p↓r$ are distinct primes, and $\mu (n)
= 0$ if $n$ is divisible by the square of a prime. Show that
$q↓n = \sum ↓{k≥1} \mu (k)\lfloor n/k\rfloor ↑2$.
|qleft \qquad c) As a consequence of (b), prove that lim$↓{n→∞}
q↓n/n↑2 = \sum ↓{k≥1}\mu (k)/k↑2$.
|qleft \qquad d) Prove that $(\sum ↓{k≥1}\mu (k)/k↑2)(\sum ↓{m≥1}/m↑2)
= 1.\xskip {\sl Hint:} When the series are absolutely convergent
we have
$$$\left(\sum ↓{k≥1} a↓k/k↑2}\left(\sum ↓{m≥1} b↓m/m↑2} = \sum
↓{n≥1} \left(\sum ↓{d\rslash n} a↓db↓{n/d}}\left/n↑2$.

\exno 11. [M22] What is the probability
that gcd($u, v) ≤ 3?$ (See Theorem C\null.) What is the {\sl average}
value of gcd($u, v)?$

\exno 12. [M24] (E. Ces\`aro.)\xskip If $u$ and $v$ are random
positive integers, what is the average number of (positive)
divisors they have in common?\xskip [{\sl Hint:} See the identity
in exercise 10(d), with $a↓k = b↓m = 1.]$
|qleft ???|newtype 58320---Compute
%folio 439 galley 10 WARNING: Some bad spots. (C) Addison-Wesley 1978	-
r Programming\quad (A.-W./Knuth)\quad ch. 4\quad f. 439\quad
g. 10
\exno 13. [HM23] Given that $u$
and $v$ are random {\sl odd} positive integers, show that they
are relatively prime with probability $8/π↑2$.

\exno 14. [M26] What are the values of $v↓1$ and $v↓2$ when
Algorithm X terminates?

\trexno 15. [M25] Design an algorithm to {\sl divide u by v modulo
m}, given positive integers $u, v$, and $m$, with $v$ relatively
prime to $m$. In other words, your algorithm should find $w$,
in the range $0 ≤ w < m$, such that $u ≡ vw$ (modulo $m)$.

\exno 16. [21] Use the text's method to find a general solution
in integers to the following sets of equations:
$$\qquad a) |tab 3$x + 7y + 11z = 1\qquad \qquad \qquad b) |tab
3x + 7y + 11z = 1⊗\3skip ⊗5x + 7y - 5z = 3⊗5x + 7y - 5z = -3\cr

\trexno 17. [M24] Show how Algorithm L
can be extended (as Algorithm A was extended to Algorithm X)
to obtain solutions of (15) when $u$ and $v$ are large.

\exno 18. [M37] Let $u$ and $v$ be odd integers, independently
and uniformly distributed in the ranges $2↑m < u < 2↑{m+1},
2↑n < v < 2↑{n+1}$. What is the {\sl exact} probability that
a single ``subtract and shift'' cycle in Algorithm B\null, namely,
an operation that starts at step B6 and then stops after step
B5 is finished, reduces $u$ and $v$ to the ranges $2↑{m↑} <
u < 2↑{m↑+1}, 2↑{n↑} < v < 2↑{n↑+1}$, as a function of $m, n,
m↑\prime $, and $n↑\prime ?$ (This exercise gives more accurate
values for the transition probabilities than in the text's model.)

\exno 19. [M24] Complete the text's derivation of (38) by establishing
(37).

\exno 20. [M26] Let $λ = \lfloor$ lg gcd($u, v)\rfloor $. Show
that the lattice-point model gives $λ = 1$ with probability
${1\over 5}$, $λ = 2$ with probability ${1\over 10}$, $λ = 3$
with probability ${1\over 20}$, etc., plus correction terms
that go rapidly to zero as $u, v$ approach infinity; hence
the average value of $λ$ is approximately ${4\over 5}$.\xskip
[{\sl Hint:} Consider the relation between the probability of
a path from $(m, n)$ to $(k + 1, k + 1)$ and a corresponding
path from $(m - k, n - k)$ to $(1, 1)$.]

\exno 21. [M25] Let $C↓{mn}, D↓{mn}$ be
the average number of subtraction and shift cycles respectively,
in Algorithm B when $u$ and $v$ are odd, \lfloor lg $u\rfloor
= m, \lfloor$ lg $v\rfloor = n$. Show that for fixed $n, C↓{mn}
= {1\over 2}m + O(1)$ and $D↓{mn} = m + O(1)$ as $m → ∞$.

\exno 22. [23] Prove Eq.\ (44).

\exno 23. [M28] Show that if $C↓{mn} = αm + βn + \gamma$ for
some constants $α, β$, and $\gamma $, then
$$$\sum ↓{1≤n≤m≤N}(N - m)(N - n)2↑{m+n-2}C↓{mn} |tab = 2↑{2N}|newtype
(|newtype {11\over 27}(α + β)N + O(1)\bigrp ,⊗\6skip \hfill
\sum ↓{1≤n≤N}(N - n)↑22↑{2n-2}C↓{nn} ⊗= 2↑{2N}\biglp {5\over
27}(α + β)N + O(1)}h12})|newtype .\cr

\trexno 24. [M30] If $v = 1$ during
Algorithm B\null, while $u$ is large, it may take fairly long for
the algorithm to determine that gcd($u, v) = 1$. Perhaps it
would be worth while to add a test at the beginning of step
B5: ``If $t = 1$, the algorithm terminates with 2$↑k$ as the
answer.'' Explore the question of whether \3skip {

\trexno 25. [M25] (R. P. Brent.)\xskip Let $u↓n, v↓n$ be the values of
$u$ and $v$ after $n$ iterations of steps B3--B5; let $X↓n =
u↓n/v↓n$, and assume that $F↓n(x)$ is the probability that $X↓n
≤ x$, for $0 ≤ x < ∞$.\xskip (a) Express $F↓{n+1}(x)$ in terms of
$F↓n(x)$, under the assumption that step B4 always branches
to B3 with probability ${1\over 2}$.\xskip (b) Let $G↓n(x) = F↓n(x)
+ 1 - F↓n(x↑{-1})$ be the probability that $Y↓n ≤ x$, for $0
≤ x ≤ 1$, where $Y↓n =$ min($u↓n, v↓n)/$max($u↓n, v↓n)$. Express
$G↓{n+1}$ in terms of $G↓n$.\xskip (c) Express $H↓n(x) =$ Pr(max($u↓{n+1},
v↓{n+1})/$max($u↓n, v↓n) < x)$ in terms of $G↓n$.

\exno 26. [M23] What is the length of the longest path from
$(m, n)$ to (0, 0) in the lattice-point model?

\trexno 27. [M28] Find values of $u, v$ with \lfloor lg $u\rfloor
= m, \lfloor$ lg $v\rfloor = n$, such that Algorithm B requires
$m + 1$ subtraction steps given $m ≥ n ≥ 1$.

|qleft \|≡3|≡3|≡.|9|4|ε[M|*/|↔M|↔o|\] |πAnalyze V. C. Harris's ``binary Euclidean algorithm.''|'{A3skip

\trexno 34. [M32] (R. W. Gosper.)\xskip Demonstrate how
to modify Algorithm B for large numbers using ideas analogous
to those in Algorithm L.

\trexno 35. [M28] (V. R. Pratt.)\xskip Extend Algorithm B to an Algorithm Y that is
analogous to Algorithm X.

\exno 36. [HM49] Find a rigorous proof that Brent's model describes the asymptotic
behavior of Algorithm B.

%folio 441 galley 11 (C) Addison-Wesley 1978	-
ch. 4\quad f. 441\quad g. 11
|qleft \18skip |newtype |raiseordrop ??{\:r 4.5.3. \xskip Analysis
of Euclid's Algorithm

|qleft }\6skip The execution time of Euclid's algorithm depends
on $T$, the number of times the division step A2 is performed.
(See Algorithm 4.5.2A and Program 4.5.2A\null.) The same quantity
$T$ is an important factor in the running time of other algorithms,
for example, the evaluation of functions satisfying a reciprocity
formula (see Section 3.3.3). We shall see in this section that
the mathematical analysis of this quantity $T$ is interesting
and instructive.

\subsectionbegin{Relation to continued fractions}
Euclid's algorithm is intimately connected with {\sl continued
fractions}, which are expressions of the form

|qleft \12skip Continued fractions have a beautiful theory that
is the subject of several books, e.g., O. Perron, {\sl Die Lehre
von den Kettenbr|spose 4uchen}, 3rd ed.\ (Stuttgart: Teubner,
1954), 2 vols.; A. Khinchin, {\sl Continued Fractions}, tr.\
by Peter Wynn (Groningen: P. Noordhoff, 1963); H. S. Wall, {\sl
Analytic Theory of Continued Fractions} (New York: Van Nostrand,
1948); see also J. Tropfke, {\sl Geschichte der Elementar-Mathematik
\bf 6} (Berlin: Gruyter, 1924), 74--84, for the early history
of the subject. It is necessary to limit ourselves to a comparatively
brief treatment of the theory here, studying only those aspects
that give us more insight into the behavior of Euclid's algorithm.

The continued fractions that are of primary interest
to us are those in which all the $b$'s in (1) are equal to unity.
For convenience in notation, let us define
$$\bslash $x↓1, x↓2, \ldotss , x↓n\bslash = 1/\biglp x↓1 +
1/(x↓2 + 1/(\cdots + 1/(x↓n) \cdot \cdot \cdot ))\bigrp .\eqno
(2)$
|qctr \6skip If $n = 0$, the symbol \bslash $x↓1, \ldotss ,
x↓n\bslash$ is taken to mean 0. Thus, for example,
$$\bslash $x↓1\bslash = {1\over x↓1} ,\qquad \bslash x↓1,
x↓2\bslash = {1\over x↓1 + (1/x↓2)} = {x↓2\over x↓1x↓2 + 1}
.\eqno (3)$
|qctr \6skip Let us also define the polynomials $Q↓n(x↓1, x↓2,
\ldotss , x↓n)$ of $n$ variables, for $n ≥ 0$, by the rule
$$$$
$$
|qleft ⊗1,\eqno if\quad $n = 0;\cr
\4skip Q↓n(x↓1, x↓2, \ldotss , x↓n) = \quad x↓1,\eqno$ if\quad
$n = 1$;

|qleft \4skip ⊗x$↓1Q↓{n-1}(x↓2, \ldotss , x↓n) + Q↓{n-2}(x↓3,
\ldotss , x↓n),\eqno$ if\quad $n > 1.\cr
\3skip (4)$
|qright \6skip Thus $Q↓2(x↓1, x↓2) = x↓1x↓2 + 1, Q↓3(x↓1, x↓2,
x↓3) = x↓1x↓2x↓3 + x↓1 + x↓3$, etc. In general, as noted by
L. Euler in the eighteenth century, $Q↓n(x↓1, x↓2, \ldotss ,
x↓n)$ is the sum of all terms obtainable by starting with $x↓1x↓2
\ldotsm x↓n$ and deleting zero or more nonoverlapping pairs of
consecutive variables $x↓jx↓{j+1}$; there are $F↓{n+1}$ such
terms. The polynomials defined in (4) are called ``continuants.``

The basic property of the $Q$-polynomials is that
$$\bslash $x↓1, x↓2, \ldotss , x↓n\bslash = Q↓{n-1}(x↓2, \ldotss
, x↓n)/Q↓n(x↓1, x↓2, \ldotss , x↓n),\qquad n ≥ 1.\eqno (5)$
|qctr \6skip This can be proved by induction, since it implies
that
$$$x↓0 + \bslash x↓1, \ldotss , x↓n\bslash = Q↓{n+1}(x↓0, x↓1,
\ldotss , x↓n)/Q↓n(x↓1, \ldotss , x↓n)$;
|qctr \6skip hence $\bslash x↓0, x↓1, \ldotss , x↓n\bslash$
is the reciprocal of the latter quantity.

The $Q$-polynomials are symmetrical in the sense
that
$$$Q↓n(x↓1, x↓2, \ldotss , x↓n) = Q↓n(x↓n, \ldotss , x↓2, x↓1).\eqno
(6)$
|qctr \6skip This follows from Euler's observation above, and
as a consequence we have
$$$Q↓n(x↓1, \ldotss , x↓n) = x↓nQ↓{n-1}(x↓1, \ldotss , x↓{n-1})
+ Q↓{n-2}(x↓1, \ldotss , x↓{n-2}).\eqno (7)$
|qctr \6skip The $Q$-polynomials also satisfy the important
identity
$$$\quad Q↓N(x↓1, \ldotss , x↓n)Q↓n(x↓2, \ldotss , x↓{n+1}) -
Q↓{n+1}(x↓1, \ldotss , x↓{n+1})Q↓{n-1}(x↓2, \ldotss , x↓n)$

|qleft \4skip = (-1)$↑n,\qquad n ≥ 1.\quad (8)$
|qright \6skip (See exercise 4.) The latter equation in connection
with (5) implies that
$$\quad $\bslash x↓1, \ldotss , x↓n\bslash = {1\over q↓0q↓1}
- {1\over q↓1q↓2} + {1\over q↓2q↓3} - \cdots + {(-1)↑{n-1}\over
q↓{n-1}q↓n} $,
|qleft where\qquad $q↓k = Q↓k(x↓1, \ldotss , x↓k).\quad (9)$
|qright \6skip Thus the $Q$-polynomials are intimately related
to continued fractions.

Every real number $X$ in the range $0 ≤ X < 1$
has a {\sl regular continued fraction} defined as follows: Let
$X↓0 = X$, and for all $n ≥ 0$ such that $X↓n ≠ 0$ let
$$$A↓{n+1} = \lfloor 1/X↓n\rfloor ,\qquad X↓{n+1} = 1/X↓n -
A↓{n+1}.\eqno (10)$
|qctr \6skip If $X↓n = 0$, the quantities $A↓{n+1}$ and $X↓{n+1}$
are not defined, and the regular continued fraction for $X$
is $\bslash A↓1, \ldotss , A↓n\bslash $. If $X ≠ 0$, this definition
guarantees that $0 ≤ X↓{n+1} < 1$, so each of the $A$'s is a
positive integer. The definition (10) clearly implies that
|qleft $\6skip X = X↓0 = {1\over A↓1 + X↓1} = {1\over A↓1 +
1/(A↓2 + X↓2)} = \cdotss$;
|qctr \6skip hence
$$$X = \bslash A↓1, \ldotss , A↓{n-1}, A↓n + X↓n\bslash \eqno
(11)$
|qctr \6skip for all $n ≥ 1$, whenever $X↓n$ is defined. Therefore,
if $X↓n = 0, X = \bslash A↓1, \ldotss , A↓n\bslash $. If $X↓n
≠ 0, X$ lies {\sl between \bslash A↓1, \ldotss , A↓n\bslash}
and $\bslash A↓1, \ldotss , A↓n + 1\bslash $, since by (7)
the quantity $q↓n = Q↓n(A↓1, \ldotss , A↓n + X↓n)$ increases
monotonically from $Q↓n(A↓1, \ldotss , A↓n)$ up to $Q↓N(A↓1,
\ldotss , A↓n + 1)$ as $X↓n$ increases from 0 to 1, and by (9)
the continued fraction increases or decreases when $q↓n$ increases,
according as $n$ is even or odd.. In fact,
$$|newtype |$X - \bslash A↓1, \ldotss , A↓n\bslash | |tab =
|\bslash A↓1, \ldotss , A↓n + X↓n\bslash - \bslash A↓1, \ldotss
, A↓n\bslash |⊗\3skip ⊗= |\bslash A↓1, \ldotss , A↓n, 1/X↓n\bslash
- \bslash A↓1, \ldotss , A↓n\bslash |\cr
\3skip ⊗= |bigabs {Q↓n(A↓2, \ldotss , A↓n, 1/X↓n)\over Q↓{n+1}(A↓1,
\ldotss , A↓n, 1/X↓n)} - {Q↓{n-1}(A↓2, \ldotss , A↓n\over Q↓n(A↓1,
\ldotss , A↓n)}|bigabs \cr
\3skip ⊗= 1/Q↓n(A↓1, \ldotss , A↓n)Q↓{n+1}(A↓1, \ldotss , A↓n,
1/X↓n)\cr
\3skip ⊗≤ 1/Q↓n(A↓1, \ldotss , A↓n)Q↓{n+1}(A↓1, \ldotss , A↓n,
A↓{n+1})\eqno (12)\cr
\6skip$ by (5), (8), and (10). Therefore \bslash $A↓1, \ldotss
, A↓n\bslash$ is an extremely close approximation to $X$. If
$X$ is irrational, it is impossible to have $X↓n = 0$ for any
$n$, so the regular continued fraction expansion in this case
is an $infinite continued fraction \bslash A↓1, A↓2, A↓3, .
. .\bslash $. The value of such a continued fraction is defined
to be lim$↓{n→∞}\bslash A↓1, A↓2, \ldotss , A↓n\bslash $, and
from the inequality (12) it is clear that this limit equals
$X$.

The regular continued fraction expansion of real
numbers has several properties analogous to the representation
of numbers in the decimal system. If we use the formulas above
to compute the regular continued fraction expansions of some
familiar real numbers, we find, for example, that
$$|newtype
 |qleft \hfill ${8\over 29}$ ⊗= \bslash 3, 1, 1, 1, 2\bslash
;\cr
\-3skip \hfill ??U19}${8\over 29}$?? ⊗= \bslash 1, 1, 9, 2,
2, 3, 2, 2, 9, 1, 2, 1, 9, 2, 2, 3, 2, 2, 9, 1, 2, 1, 9, 2,
2, 3, 2, 2, 9, 1, . .3.\bslash ;

|qleft \2skip \hfill |spose $↑3\sqrt{2} ⊗= 1 + \bslash 3, 1,
5, 1, 1, 4, 1, 1, 8, 1, 14, 1, 10, 2, 1, 4, 12, 2, 3, 2, 1,
3, 4, 1, 1, 2, 14, 3, . . .\bslash ;\cr
\2skip \hfill π ⊗= 3 + \bslash 7, 15, 1, 292, 1, 1, 1, 2, 1,
3, 1, 1, 14, 2, 1, 1, 2, 2, 2, 2, 1, 84, 2, 2, 1, 1, 15, 3,
13, . . .\bslash ;\cr
\2skip \hfill e ⊗= 2 + \bslash 1, 2, 1, 1, 4, 1, 1, 6, 1, 1,
8, 1, 1, 12, 1, 1, 12, 1, 1, 14, 1, 1, 16, 1, 1, 18, 1, . .
.\bslash ;\cr
\hfill \gamma ⊗= \bslash 1, 1, 2, 1, 2, 1, 4, 3, 13, 5, 1,
1, 8, 1, 2, 4, 1, 1, 40, 1, 11, 3, 7, 1, 7, 1, 1, 5, . . .\bslash
;\cr
\2skip \hfill \phi ⊗= 1 + \bslash 1, 1, 1, 1, 1, 1, 1, 1, 1,
. . .\bslash .\eqno (13)\cr
\6skip |newtype$ The numbers $A↓1, A↓2, . . $. are called the
{\sl partial quotients} of $X$. Note the regular pattern that
appears in the partial quotients for \sqrt{8/29}, $\phi $, and
$e$; the reasons for this behavior are discussed in exercises
12 and 16. There is no apparent pattern in the partial quotients
for |spose $↑3\sqrt{2}, π$, or $\gamma $.

It is interesting to note that the ancient Greek's
first definition of real numbers, once they had discovered the
existence of irrationals, was essentially in terms of infinite
continued fractions. (Later they adopted the suggestion of Eudoxus
that $x = y$ should be defined instead as ``$x < r$ if and only
if $y < r$, for all rational $r$.'' See O. Becker, {\sl Quellen
und Studien zur Geschichte Math., Astron., Physik} (B) {\bf 2}
(1933), 311--333.)
|qleft
%folio 444 galley 12 (C) Addison-Wesley 1978	-
\def\bslash{\char'477 } \def\vbslash{\char'477017 } % boldface slashes (vol. 2 only)
\yskip When $X$ is a rational number, the regular
continued fraction corresponds in a natural way to Euclid's
algorithm. Let us assume that $X = v/u$, where $u > v ≥ 0$.
The regular continued fraction process starts with $X↓0 = X$;
let us define $U↓0 = u$, $V↓0 = v$. Assuming that $X↓n = V↓n/U↓n
≠ 0$, (10) becomes
$$\baselineskip15pt\cpile{A↓{n+1} = \lfloor U↓n/V↓n\rfloor,\cr
X↓{n+1} = U↓n/V↓n - A↓{n+1} = (U↓n\mod V↓n)/V↓n.\cr}\eqno(14)$$
Therefore, if we define
$$U↓{n+1} = V↓n,\qquad V↓{n+1} = U↓n\mod V↓n,\eqno (15)$$
the condition $X↓n = V↓n/U↓n$ holds throughout
the process. Furthermore, (15) is precisely the transformation
made on the variables $u$ and $v$ in Euclid's algorithm (see Algorithm
4.5.2A\null, step A2). For example, since ${8\over 29} = \bslash
3, 1, 1, 2\bslash$, we know that Euclid's algorithm applied
to $u = 29$ and $v = 8$ will require exactly five division steps,
and the quotients $\lfloor u/v\rfloor$ in step A2 will be successively
3, 1, 1, 1, and 2. Note that the last partial quotient $A↓n$
must be 2 or more when $X↓n = 0$, $n ≥ 1$, since $X↓{n-1}$ is
less than unity.

From this correspondence with Euclid's algorithm
we can see that the regular continued fraction for $X$ terminates
at some step with $X↓n = 0$ if and only if $X$ is rational;
for it is obvious that $X↓n$ cannot be zero if $X$ is irrational,
and, conversely, we know that Euclid's algorithm always terminates.
If the partial quotients obtained during Euclid's algorithm
are $A↓1$, $A↓2$, $\ldotss$, $A↓n$, then we have, by (5),
$${v\over u} = {Q↓{n-1}(A↓2, \ldotss , A↓n)\over Q↓n(A↓1, A↓2,
\ldotss , A↓n)} .\eqno (16)$$
This formula holds also if Euclid's algorithm is
applied for $u < v$, when $A↓1 = 0$. Furthermore, because of
(8), $Q↓{n-1}(A↓2, \ldotss , A↓n)$ and $Q↓n(A↓1, A↓2, \ldotss
, A↓n)$ are relatively prime, and the fraction on the right-hand
side of (16) is in lowest terms; therefore
$$u = Q↓n(A↓1, A↓2, \ldotss , A↓n)d,\qquad v = Q↓{n-1}(A↓2,
\ldotss , A↓n)d,\eqno (17)$$
where $d = \gcd(u, v)$.

\subsectionbegin{The worst case} We can now apply
these observations to determine the behavior of Euclid's algorithm
in the ``worst case,'' or in other words to give an upper bound
on the number of division steps. The worst case occurs when
the inputs are consecutive Fibonacci numbers:

\algbegin Theorem F ({\rm G. Lam\'e, 1845}). {\sl For $n ≥ 1$, let
$u$ and $v$ be integers with $u > v > 0$ such that Euclid's algorithm
applied to $u$ and $v$ requires exactly $n$ division steps, and such
that $u$ is as small as possible satisfying these conditions.
Then $u = F↓{n+2}$ and $v = F↓{n+1}$.}

\proofbegin By (17), we must have $u = Q↓n(A↓1,
A↓2, \ldotss , A↓n)d$ where $A↓1$, $A↓2$, $\ldotss$, $A↓n$, $d$ are positive
integers and $A↓n ≥ 2$. Since $Q↓n$ is a polynomial with nonnegative
coefficients, involving all of the variables, the minimum value
is achieved only when $A↓1=1$, $\ldotss$, $A↓{n-1}=1$, $A↓n=2$, $d=1$. Putting
these values in (17) yields the desired result.\quad\blackslug

\yyskip (This theorem has the historical claim
of being the first practical application of the Fibonacci sequence;
since then many other applications of the Fibonacci numbers
to algorithms and to the study of algorithms have been discovered.)

As a consequence of Theorem F we have an important
corollary:

\thbegin Corollary. {\sl If $0 ≤ u, v < N$, the number of division
steps required when Algorithm 4.5.2A is applied to $u$ and
$v$ is at most $\lceil log↓|≤f (\sqrt{5}N)\rceil - 2.}

\proofbegin By Theorem F\null, the maximum number
of steps, $n$, occurs when $u = F↓n$ and $v = F↓{n+1}$, where
$n$ is as large as possible with $F↓{n+1} < N$. (The first division
step in this case merely interchanges $u$ and $v$ when $n >
1.)$ Since $F↓{n+1} < N$, we have $\phi ↑{n+1}/\sqrt{5} < N
($see Eq.\ 1.2.8--15), so $n + 1 <$ log$↓|≤f (\sqrt{5}N)$.

|qleft \12skip Note that log$↓|≤f (\sqrt{5}N)$ is approximately
2.078 ln $N + 1.672 \approx 4.785$ log$↓{10} N + 1.672$. See
exercises 31 and 36 for extensions of Theorem F.

\subsectionbegin{An approximate model} Now that we
know the maximum number of division steps that can occur, let
us attempt to find the {\sl average} number. Let $T(m, n)$ be
the number of division steps that occur when $u = m, v = n$
are input to Euclid's algorithm. Thus
$$$T(m, 0) = 0;\qquad T(m, n) = 1 + T(n, m$ mod $n)\qquad$ if\qquad
$n ≥ 1.\eqno (18)$
|qctr \6skip Let $T↓n$ be the average number of division steps
when $v = n$ and when $u$ is chosen at random; since only the
value of $u$ mod $v$ affects the algorithm after the first division
step, we may write
$$$T↓n = {1\over n} \sum ↓{0≤k<n} T(k, n).\eqno (19)$
|qctr \6skip For example, $T(0, 5) = 1, T(1, 5) = 2, T(2, 5)
= 3, T(3, 5) = 4, T(4, 5) = 3$, so
$$$T↓5 = {1\over 5}(1 + 2 + 3 + 4 + 3) = 2{3\over 5}$.
|qctr \6skip \quad In order to estimate $T↓n$ for large $n$,
let us first try an approximation suggested by R. W. Floyd:
We might assume that, for $0 ≤ k < n, n$ is essentially ``random''
modulo $k$, so that we can set
$$$T↓n \approx 1 + {1\over n} (T↓0 + T↓1 + \cdots + T↓{n-1})$.
|qctr \6skip Then $T↓n \approx S↓n$, where the sequence $\langle
S↓n\rangle$ is the solution to the recurrence relation
$$$S↓0 = 0,\qquad S↓n = 1 + {1\over n} (S↓0 + S↓1 +\cdot \cdot
\cdot + S↓{n-1}),\qquad n ≥ 1.\eqno (20)$
|qctr \6skip (This approximation is analogous to the ``lattice-point
model'' used to investigate Algorithm B in Section 4.5.2.)

The recurrence (20) is readily solved by the use
of generating functions. A more direct way to solve it, analogous
to ur solution of the lattice-point model, is by noting that
$$
|qctr \hfill S$↓{n+1} ⊗= 1 + {1\over n + 1} (S↓0 + S↓1 +\cdot
\cdot \cdot + S↓{n-1} + S↓n)\cr
\6skip ⊗= 1 + {1\over n + 1} \biglp n(S↓n - 1) + S↓n|newtype
) |newtype = S↓n + {1\over n + 1} ;\cr
\6skip$ hence $S↓n$ is 1 + ${1\over 2}$ +\cdot \cdot \cdot +
1/$n = H↓n$, a harmonic number. Therefore the approximation
$T↓n \approx S↓n$ would tell us that $T↓n \approx$ ln $n + O(1)$.

Comparison of this approximation with tables of
the true value of $T↓n$ show, however, that ln $n$ is too large;
$T↓n$ does not grow this fast. One way to account for the fact
that this approximation is too pessimistic is to observe that
the average value of $n$ mod $k$ is less than the average value
of ${1\over 2}k$, in the range $1 ≤ k ≤ n:$
$$${1\over n}$ \sum ↓{1≤k≤n} (n mod $k) |tab = {1\over n} \sum
??1≤q≤n??\lfloor n/(q+1)\rfloor <k≤\lfloor n/q\rfloor ?? (n
- qk) = n - {1\over n} \sum ↓{1≤q≤n} q\left(\left({\lfloor n/q\rfloor
+ 1\atop 2}} - \left({\lfloor n/(q + 1)\rfloor + 1\atop 2}}}⊗\6skip
⊗= n - {1\over n} \sum ↓{1≤q≤n} \left({\lfloor n/q\rfloor +
1\atop 2}} = \left(1 - {π↑2\over 12}}n + O($log $n)\eqno (21)\cr
\6skip $\biglp$cf.\ exercise 4.5.2--10(c)$\bigrp$. This is only
about .1775n$, not $.25n$, so the value of $n$ mod $k$ tends
to be smaller than the above model predicts; Euclid's algorithm
works faster than we might expect.
|qleft ??????|newtype 58320
%folio 448 galley 13 (C) Addison-Wesley 1978	-
\def\bslash{\char'477 } \def\vbslash{\char'477017 } % boldface slashes (vol. 2 only)
\subsectionbegin{A continuous model} The behavior
of Euclid's algorithm with $v = N$ is essentially determined
by the behavior of the regular continued fraction process when
$X = 0/N$, $1/N$, $\ldotss$, $(N - 1)/N$. Assuming that $N$ is very
large, we are led naturally to a study of regular continued
fractions when $X$ is a random real number uniformly distributed in
$[\,0, 1)$. Therefore let us define the distribution function
$$F↓n(x) =\hjust{probability that }X↓n ≤ x,\qquad\hjust{for }0
≤ x ≤ 1,\eqno (22)$$
given a uniform distribution of $X = X↓0$. By the
definition of regular continued fractions, we have $F↓0(x) =
x$, and
$$\eqalignno{F↓{n+1}(x) ⊗= \sum ↓{k≥1}\hjust{probability that}(k
≤ 1/X↓n ≤ k + x)\cr
⊗= \sum ↓{k≥1}\hjust{probability that}\biglp 1/(k + x) ≤ X↓n ≤
1/k\bigrp\cr
⊗= \sum ↓{k≥1} \biglp F↓n(1/k) - F↓n\biglp 1/(k +x)\bigrp\bigrp.⊗(23)\cr}$$
If the distributions $F↓0(x)$, $F↓1(x)$, $\ldots$ defined
by these formulas approach a limiting distribution $F↓∞(x) =
F(x)$, we will have
$$\chop to 12pt{F(x) = \sum ↓{k≥1} \biglp F(1/k) - F\biglp 1/(k + x)\bigrp\bigrp.}
\eqno (24)$$
One function that satisfies this relation
is $F(x) = \log↓b (1 + x)$, for any base $b > 1$; see exercise
19. The further condition $F(1) = 1$ implies that we should
take $b = 2$. Thus it is reasonable to make a guess that $F(x)
= \lg (1 + x)$, and that $F↓n(x)$ approaches this behavior.

We might conjecture, for example, that $F({1\over
2}) =\lg({3\over 2}) \approx 0.58496$; let us see how close
$F↓n({1\over 2})$ comes to this value for small $n$. We have
$F↓0({1\over 2})={1\over2}$, and
|qleft $$
|qctr $\hfill F↓1({1\over 2}) ⊗= {1\over 1} - {1\over 1 + {1\over
2}?? + {1\over 2} - {1\over 2 + {1\over 2}?? +\cdot \cdot \cdot \cr
\5skip ⊗= 2\left({1\over 2} - {1\over 3} + {1\over 4} - {1\over
5} +\cdots } = 2(1 - ln 2) \approx 0.6137;\cr
$\5skip \hfill F↓2({1\over 2}) ⊗= \sum ↓{m≥1} {2\over m} \left({1\over
2m + 2} - {1\over 3m + 2} + {1\over 4m + 2} - {1\over 5m + 2}
+\cdot \cdot \cdot }\cr
\4skip ⊗= \sum ↓{m≥1} {2\over m↑2} \left({1\over 2} - {1\over
3} + {1\over 4} -\cdot \cdot \cdot }\cr
\4skip ⊗ - \sum ↓{m≥1} {4\over m} \left({1\over 2m(2m + 2)}
- {1\over 3m(3m + 2) +\cdot \cdot \cdot }\cr
\5skip ⊗= {1\over 3}π↑2(1 - ln 2) - \sum ↓{m≥1} ${4S↓m\over
m↑2}$ ,\cr
\6skip for some positive constant $A$. The error term was improved
to $O(e↑{-An})$ by Paul L|spose 1evy shortly afterward [{\sl
Bull.\ Soc.\ Math.\ de France \bf 57} (1929), 178--194];?? but
Gauss's problem, namely to find the asymptotic behavior of $F↓n(x)
- lg(1 + $x)$, was not really resolved until 1974, when Eduard
Wirsing published a beautiful analysis of the situation [{\sl
Acta Arithmetica \bf 24} (1974), 507--528]. We shall study the
simplest aspects of Wirsing's approach here, since his method
appears to be the best way to handle problems of this type.

If $G$ is any function of $x$ defined for $0≤ x
≤ 1$, let $SG$ be the function defined by
$$$SG(x) = \sum ↓{k≥1} \left(G\left({1\over k}} - G\left({1\over
k + x}}}.\eqno (25)$
|qctr \6skip Thus, $S$ is an operator that changes one function
into another. In particular, by (23) we have $F↓{n+1}(x) = SF↓n(x)$,
hence
$$$F↓n = S↑nF↓0.\eqno (26)$
|qctr \6skip (In this discussion $F↓n$ stands for a distribution
function, {\sl not} for a Fibonacci number.) Note that $S$ is
a ``linear operator''; i.e., $S(cG) = c(SG)$ for constants $c$,
and $S(G↓1 + G↓2) = SG↓1 + SG↓2$.

Now if $G$ has a bounded first derivative, we can
differentiate (25) term by term to show that
$$$(SG)↑\prime (x) = \sum ↓{k≥1} {1\over (k + x)↑2} G↑\prime
\left({1\over k + x}},\eqno (27)$
|qctr \6skip ---------------

|qleft \3skip |newtype ?? An exposition of L|spose 1evy's interesting
proof appeared in the first edition of this book.

|qleft \6skip |newtype hence $SG$ also has a bounded first derivative.
(Term-by-term differentiation of a convergent series is justified
when the series of derivatives is uniformly convergent; cf.\
K. Knopp, {\sl Theory and Application of Infinite series (}Glasgow:
Blackie, 1951), 47.)

Let $H = SG$, and let $g(x) = (1 + x)G↑\prime (x),
h(x) = (1 + x)H↑\prime (x)$. It follows that
$$$h(x) = \sum ↓{k≥1} {1 + x\over (k + x)↑2} \left(1 + {1\over
k + x}↑{-1}g\left({1\over k + x} =\sum ↓{k≥1}\left({k\over k
+ 1 + x) - {k - 1\over k + x}}g\left({1\over k + x}}}$.
|qctr \6skip In other words, $h = Tg$, where $T$ is the linear
operator defined by
$$$Tg(x) = \sum ↓{k≥1} \left({k\over k + 1 + x} - {k - 1\over
k + x}}g\left({1\over k + x}}.\eqno (28)$
|qctr \quad Continuing, we see that if $g$ has a bounded first
derivative, we can differentiate term by term to show that $Tg$
does also:
$$(Tg)↑\prime (x) |tab = - \sum ↓{k≥1} \left(${k\over (k + 1
+ x)↑2}$ \left(g\left(${1\over k + x}$} - g\left(${1\over k
+ 1 + x}$}} + ${1 + x\over (k + x)↑3(k + 1 + x)}$ g↑\prime \left(${1\over
k + x}$}}.⊗\6skip ⊗= - \sum ↓{k≥1} \left(${k\over (k + 1 + x)↑2}$
\left(g\left(${1\over k + x}$} - g\left(${1\over k + 1 + x}$}}
+ ${1 + x\over (k + x)↑3(k + 1 + x)}$ g↑\prime \left(${1\over
k + x}$}}.\cr
\6skip There is consequently a third linear operator, $U$, such
that $(Tg)↑\prime = -U(g)↑\prime )$, namely
$$$U\varphi (x) = \sum ↓{k≥1} \left({k\over (k + 1 + x)↑2} \int
↑{1/(k+x)}↓{1/(k+1+x)} \varphi (t)dt + {1 + x\over (k + x)↑3(k
+ 1 + x)} \varphi \left({1\over k + x}}}$.
|qctr (29)
|qright \6skip ????????????|newtype 58320---Computer Programming\quad
%folio 451 galley 14 (C) Addison-Wesley 1978	-
\def\bslash{\char'477 } \def\vbslash{\char'477017 } % boldface slashes (vol. 2 only)
What is the relevance of all this to our problem?
Well, if we set
$$\eqalignno{F↓n(x)⊗=\lg(1 + x) + R↓n\biglp\lg(1 + x)\bigrp,⊗(30)\cr
\noalign{\vskip4pt}
f↓n(x) = (1 + x)F↑\prime↓{\!n}(x) = {1\over
\ln 2} \biglp 1 + R↑\prime↓{\!n}\biglp\lg(1 + x)\bigrp\bigrp,⊗(31)\cr}$$
,\eqno (31)\cr
we have
$$f↑\prime↓{\!n}(x) = R↓{\!n}↑{\prime\prime}\biglp\lg(1 + x)\bigrp/\biglp(\ln 2)↑2(1
+ x)\bigrp;\eqno (32)\cr$$
the effect of the $\lg(1 + x)$ term disappears, after
these transformations. Furthermore since $F↓n = S↑nF↓0$ we have
$f↓n = T↑nf↓0$ and $f↑\prime↓{\!n} = (-1)↑nU↑nf↑\prime↓{0}$,
and $F↓n$ and $f↓n$ have bounded derivatives, by induction on
$n$. Thus (32) becomes
$$(-1)↑nR↓{\!n}↑{\prime\prime}\biglp\lg(1+x)\bigrp=(1+x)(\ln2)↑2U↑nf↓0↑\prime(x).
\eqno (33)$$
Now $F↓0(x) = x$, $f↓0(x) = 1 + x$, and
$f↓0↑\prime(x)$ is the constant function 1. We are going
to show that the operator $U↑n$ takes the constant function
into a function with very small values, hence $|R↓{\!n}↑{\prime\prime}(x)|$ must
be very small for $0 ≤ x ≤ 1$. Finally we can clinch the argument
by showing that $R↓n(x)$ itself is small: Since $R↓n(0) = R↓n(1)
= 1$, it follows from a well-known interpolation formula (cf.\
exercise 4.6.4--15 with $x↓0 = 0$, $x↓1 = x$, $x↓2 = 1$) that
$$R↓n(x) = -{x(1 - x)\over 2} R↓{\!n}↑{\prime\prime}\biglp\xi(x)\bigrp\eqno (34)$$
for some function $\xi (x)$, where $0 ≤ \xi (x) ≤ 1$ when
$0 ≤ x ≤ 1$.

Thus everything hinges on our being able to prove
that $U↑n$ produces small function values, where $U$ is the
linear operator defined in (29). Note that $U$ is a {\sl positive}
operator in the sense that $U\varphi (x) ≥ 0$ for all $x$ if
$\varphi (x) ≥ 0$ for all $x$. It follows that $U$ is order-preserving:
If $\varphi ↓1(x) ≤ \varphi ↓2(x)$ for all $x$ then $U\varphi
↓1(x) ≤ U\varphi ↓2(x)$ for all $x$.

One way to exploit this property is to find a function
$\varphi$ for which we can calculate $U\varphi$ exactly and
to use constant multiples of this function to bound the ones
we are really interested in. First let's look for a function $g$
such that $Tg$ is easy to compute. If we consider functions
defined for all $x ≥ 0$ instead of only on $[\,0, 1]$, it is easy
to remove the summation from (25) by observing that
$$SG(x + 1) - SG(x) = G\left(1\over 1 + x\right) - \lim↓{k→∞}
G\left(1\over k + x\right) = G\left(1\over 1 + x\right) - G(0)\eqno(35)$$
when $G$ is continuous. Since $T\biglp (1 + x)G↑\prime
\bigrp = (1 + x)(SG)↑\prime $, it follows (see exercise 20)
that
$${Tg(x)\over x + 1} - {Tg(x + 1)\over x + 2} = \left({1\over
x + 1} - {1\over x + 2}}g\left(1\over 1 + x\right).\eqno (36)$$
If we set $Tg(x) = 1/(x + 1)$, we find that the
corresponding value of $g(x)$ is $1 + x - 1/(1 + x)$. Let $\varphi
(x) = g↑\prime (x) = 1 + 1/(1 + x)↑2$, so that $U\varphi (x)
= -(Tg)↑\prime (x) = 1/(1 + x)↑2$; this is the function $\varphi$
we have been looking for.

For this choice of $\varphi$ we have $2 ≤ \varphi
(x)/U\varphi (x) = (1 + x)↑2 + 1 ≤ 5$ for $0 ≤ x ≤ 1$, hence
$$\textstyle{1\over 5}\varphi ≤ U\varphi ≤ {1\over 2}\varphi .$$
By the positivity of $U$ and $\varphi$ we can apply
$U$ to this inequality again, obtaining ${1\over 25}\varphi
≤ {1\over 5}U\varphi ≤ U↑2\varphi ≤ {1\over2}U\varphi ≤{1\over 4}\varphi
$; and after $n - 1$ applications we have
$$5↑{-n}\varphi ≤ U↑n\varphi ≤ 2↑{-n}\varphi \eqno (37)$$
for this particular $\varphi $. Let $\psi (x) =
f↑\prime↓{0}(x) = 1$ be the constant function; then for $0
≤ x ≤ 1$ we have ${5\over 4}\psi ≤ \varphi ≤ 2\psi $, hence
$$\textstyle{5\over 8}5↑{-n}\psi ≤ {1\over 2}5↑{-n}\varphi ≤ {1\over
2}U↑n\varphi ≤ U↑n\psi ≤ {4\over 5}U↑n\varphi ≤ {4\over 5}2↑{-n}\varphi
≤ {8\over 5}2↑{-n}\psi .$$
It follows by (33) that
$${\textstyle{5\over 8}}(\ln 2)↑25↑{-n} ≤ (-1)↑nR↓{\!n}↑{\prime\prime}(x)
 ≤ {\textstyle{16\over 5}}(\ln
2)↑22↑{-n},\qquad 0 ≤ x ≤ 1,$$
hence by (30) and (34) we have proved the following result:

\thbegin Theorem W.  $F↓n(x) = \lg(1 + x) + O(2↑{-n})$.
{\sl In fact, $F↓n(x) - \lg(1 + x) lies between ${5\over 16}(-1)↑{n+1}5↑{-n}\biglp
\ln(1 + x)\bigrp\biglp\ln 2/(1 + x))$ and ${8\over 5}(-1)↑{n+1}2↑{-n}\biglp
\ln(1 + x)\bigrp\biglp\ln 2/(1 + x)\bigrp$, for $0 ≤ x ≤ 1$.\quad\blackslug

With a slightly different choice of $\varphi
$, we can obtain tighter bounds (see exercise 21). In fact,
Wirsing went much further in his paper, proving that
$$F↓n(x) = \lg(1 + x) + (-λ)↑n\Psi (x) + O\biglp x(1 - x)(λ - 0.031)↑n\bigrp,
\eqno(38)$$
where
$$\baselineskip15pt\eqalign{λ ⊗= 0.30366\ 30028\ 98732\ 65860\ \ldots\cr
⊗= \bslash 3, 3, 2, 2, 3, 13, 1, 174, 1, 1, 1, 2, 2,
2, 1, 1, 1, 2, 2, 1, \ldotss\bslash\cr} \eqno (39)$$
is a fundamental constant (apparently unrelated to more familiar constants), and
where $\Psi$ is an interesting function that is analytic in
the entire complex plane except for the negative real axis from
$-1$ to $-∞$. Wirsing's function satisfies $\Psi (0) = \Psi (1)
= 0$, $\Psi ↑\prime (0) < 0$, and $S\Psi = -λ\Psi $; thus by (35)
it satisfies the identity
$$\Psi (z) - \Psi (z + 1) = {1\over λ} \Psi \left(1\over 1+ z\right).\eqno (40)$$
Furthermore, Wirsing demonstrated that
$$\Psi \left(-{u\over v} + {i\over N}\right) = cλ↑{-n}\log N +
O(1)\qquad\hjust{as }N → ∞,\eqno (41)$$
where $c$ is a constant and $n = T(u, v)$ is the
number of iterations when Euclid's algorithm is applied to the
integers $u > v > 0$.

A complete solution to Gauss's problem was found a few years later by K. I. Babenko
[{\sl Doklady Akad.\ Nauk CCCP \bf238} (1978), 1021--1024],
who used powerful techniques of functional analysis to prove that
$$F↓n(x)=\lg(1+x)+\sum↓{j≥2}c↓jλ↓j↑n\Psi↓j(x)\eqno(42)$$
for all $0≤x≤1$, $n≥1$. Here $|λ↓2|>|λ↓3|≥|λ↓4|≥\cdots$, and each $\Psi↓j(z)$
is an analytic function in the complex plane except for a cut at
$[-∞,-1]$. The function $\Psi↓2$ is Wirsing's $\Psi$, and $λ↓2=-λ$, while
$λ↓3=????$ (to be supplied). Babenko also
established further properties of the $λ↓j$,
proving in particular that the sum for $j≥k$ in (42) is bounded by
$(π↑2/6)|λ↓k|↑{n-1}\min(x,1-x)$.
%folio 460 galley 15 (C) Addison-Wesley 1978	-
\def\bslash{\char'477 } \def\vbslash{\char'477017 } % boldface slashes (vol. 2 only)
\subsectionbegin{From continuous to discrete} We
have now derived results about the probability distributions
for continued fractions when $X$ is a real number uniformly
distributed in the interval $[\,0, 1)$. But since a real number
is rational with probability zero (almost all numbers are irrational),
these results do not apply directly to Euclid's algorithm. Before
we can apply Theorem W to our problem, some technicalities must
be overcome. Consider the following observation based on elementary
measure theory:

\thbegin Lemma M. {\sl Let $I↓1$, $I↓2$, $\ldotss$, $J↓1$, $J↓2$,
$\ldots$ be pairwise disjoint intervals contained in the interval
$[\,0, 1)$, and let $\Iscr = \union↓{k≥1} I↓k$, $\Jscr = \union↓{k≥1}
J↓k$. Assume that $\Kscr = [\,0, 1] \rslash (\Iscr ∪ \Jscr)$ has measure zero.
Let $P↓n$ be the set $\{0/n, 1/n, \ldotss , (n - 1)/n\}$. Then}
$$\lim↓{n→∞} {\| \Iscr ∩ P↓n\| \over n} = \mu (\Iscr).\eqno (43)$$
Here $\mu (\Iscr)$ is the Lebesgue measure of $\Iscr$,
namely, $\sum ↓{k≥1}\hjust{length}(I↓k)$; and ${\|\Iscr ∩ P↓n\|}$
denotes the number of elements of $\Iscr ∩ P↓n$.

\proofbegin Let $\Iscr↓N = \union↓{1≤k≤N}
I↓k$, $\Jscr↓N = \union↓{1≤k≤N} J↓k$. Given $ε > 0$, find $N$
large enough so that $\mu (\Iscr↓N) + \mu (\Jscr↓N) ≥ 1 - ε$, and
let
$$\Kscr↓N = \Kscr ∪ \union↓{k>N} I↓k ∪ \union↓{k>N} J↓k.$$
If $I$ is an interval, having any of the forms
$(a, b)$ or $[a, b)$ or $(a, b]$ or $[a, b]$, it is clear that
$\mu (I) = b - a$ and
$$n\mu (I) - 1 ≤ \| I ∩ P↓n\| ≤ n\mu (I) + 1.$$
|qctr \6skip Now let $r↓n = \parallel \Iscr↓N ∩ P↓n\parallel ,
s↓n = \parallel \Jscr↓N ∩ P↓n\parallel , t↓n = \parallel \Kscr↓N
∩ P↓n\parallel $; we have
$$$r↓n + s↓n + t↓n = n$;
|qctr \6skip n\mu (\Iscr$↓N) - N ≤ r↓n ≤ n\mu (\Iscr↓N) + N$;
|qctr \3skip n\mu (\Iscr$↓N) - N ≤ s↓n ≤ n\mu (\Jscr↓N) + N$.
|qctr \6skip Hence
$$$\mu (\Iscr) - {N\over n} - ε ≤ \mu (\Iscr↓N) - {N\over n} ≤ {r↓n\over
n} ≤ {r↓n + t↓n\over n}$

|qleft \4skip = 1 - ${s↓n\over n}$ ≤ 1 - \mu (\Jscr$↓N) + {N\over
n} ≤ \mu (\Iscr) + {N\over n} + ε$.
|qright \6skip This holds for all $n$ and for all $ε$; hence
lim$↓{n→∞} r↓n/n = \mu (\Iscr)$.

|qleft \12skip \quad Exercise 25 shows that Lemma M is not trivial,
in the sense that some rather restrictive hypotheses are needed
to prove (43).

\subsectionbegin{Distribution of partial quotients}
Now we can put Theorem W and Lemma M together to derive some
solid facts about Euclid's algorithm.

\thbegin Theorem E. {\sl Let $n$ and $k$ be positive
integers, and let p↓k(a, n) be the probability that the (k +
1)st quotient A↓{k+1} in Euclid's algorithm is equal to a, when
v = n and u is chosen at random. Then}
$$\mathop{lim↓{$n→∞} p↓k(a, n) = F↓k\left({1\over a}} - F↓k\left({1\over
a + 1}}$,
|qctr \6skip where F$↓k(x) is the distribution function (21)$.

|qleft \12skip Proof. \xskip The set $\Iscr$ of all $X$ in [0,
1) for which $A↓{k+1} = a$ is a union of disjoint intervals,
and so is the set $\Jscr$ of all $X$ for which $A↓{k+1} ≠ a$.
Lemma M therefore applies, with $\Kscr$ the set of all $X$ for
which $A↓{k+1}$ is undefined. Furthermore, $F↓k(1/a) - F↓k(1/(a
+ 1)|newtype ) |newtype$ is the probability that $1/(a + 1)
< X↓k ≤ 1/a$, which is $\mu (\Iscr)$, the probability that $A↓{k+1}
= a$.

|qleft \12skip \quad As a consequence of Theorems E and W\null, we
can say that a quotient equal to $a$ occurs with the approximate
probability
$$lg(1 + 1/$a) - lg\biglp 1 + 1/($a + 1)|newtype ) |newtype
=$ lg\biglp$ (a + 1)↑2/(|newtype (a + 1)↑2 - 1)\bigrp $.
|qctr \6skip Thus
$$|tab ⊗⊗a quotient of 1 occurs about lg(${4\over 3}$)⊗= 41.504
percent of the time;\cr
\3skip ⊗a quotient of 2 occurs about lg(${9\over 8}$)⊗= 16.992
percent of the time;\cr
\3skip ⊗a quotient of 3 occurs about lg(${16\over 15}$)⊗= 9.311
percent of the time;\cr
\3skip ⊗a quotient of 4 occurs about lg(${25\over 24}$)⊗= 5.890
percent of the time.\cr
\6skip Actually, if Euclid's algorithm produces the quotients
$A↓1, A↓2, \ldotss , A↓t$, the nature of the proofs above will
guarantee this behavior only for $A↓k$ when $k$ is comparatively
small with respect to $t$; the values $A↓{t-1}, A↓{t-2}, . .
$. are not covered by this proof. But we can in fact show that
the distribution of the last quotients $A↓{t-1}, A↓{t-2}, .
. $. is essentially the same as the first.

For example, consider the regular continued fraction
expansions for the fractions whose denominator is 29:
|qleft β????????\quad ??\quad ??\quad ??\quad ??\quad ??\quad
??\quad ??\quad ??\quad ??\quad ??\quad
%folio 461 galley 16 (C) Addison-Wesley 1978	-
\def\bslash{\char'477 } \def\vbslash{\char'477017 } % boldface slashes (vol. 2 only)
?????????|newtype 58320---Computer Programming\quad (A.-W./Knuth)\quad
ch. 4\quad f. 461\quad g. 16

|qleft \12skip |tab \qquad \qquad \qquad \quad |tab \qquad \qquad
\qquad \qquad \quad |tab \qquad \qquad \qquad \qquad \qquad
|tab \qquad \qquad \qquad \quad |tab |zfa ⊗${1\over 29}$ = \bslash
29\bslash ⊗${8\over 29}$ = \bslash 3, 1, 1, 1, 2\bslash ⊗${15\over
29}$ = \bslash 1, 1, 14\bslash ⊗${22\over 29}$ = \bslash
1, 3, 7\bslash \cr
${2\over 29}$ = \bslash 14, 2\bslash ⊗${9\over 29}$ = \bslash
3, 4, 2\bslash ⊗${16\over 29}$ = \bslash 1, 1, 4, 3\bslash
⊗${23\over 29}$ = \bslash 1, 3, 1, 5\bslash \cr
${3\over 29}$ = \bslash 9, 1, 2\mp ⊗${10\over 29}$ = \bslash
2, 1, 9\mp ⊗${17\over 29}$ = \bslash 1, 1, 2, 2, 2\bslash
⊗${24\over 29}$ = \bslash 1, 4, 1, 4\bslash \cr
${4\over 29}$ = \bslash 7, 4\bslash ⊗${11\over 29}$ = \bslash
2, 1, 1, 1, 3\bslash ⊗${18\over 29}$ = \bslash 1, 1, 1, 1,
1, 3\bslash ⊗${25\over 29}$ = \bslash 1, 6, 4\bslash \cr
${5\over 29}$ = \bslash 5, 1, 4\bslash ⊗${12\over 29}$ = \mp
2, 2, 2, 2\bslash ⊗${19\over 29}$ = \bslash 1, 1, 1, 9\bslash
⊗${26\over 29}$ = \bslash 1, 8, 1, 2\bslash \cr
${6\over 29}$ = \bslash 4, 1, 5\bslash ⊗${13\over 29}$ = \bslash
2, 4, 3\bslash ⊗${20\over 29}$ = \bslash 1, 2, 4, 2\bslash
⊗${27\over 29}$ = \bslash 1, 13, 2\bslash \cr
${7\over 29}$ = \bslash 4, 7\bslash ⊗${14\over 29}$ = \bslash
2, 14\bslash ⊗${21\over 29}$ = \bslash 1, 2, 1, 1, 1, 2\bslash
⊗${28\over 29}$ = \bslash 1, 28\bslash \cr
\6skip Several things can be observed in this table.

a) As mentioned earlier, the last quotient is always
2 or more. Furthermore, we have the obvious identity
$$$\bslash x↓1, \ldotss , x↓{n-1}, x↓n + 1\bslash = \bslash
x↓1, \ldotss , x↓{n-1}, x↓n, 1\bslash \eqno (44)$
|qctr \6skip which shows how partial fractions whose last quotient
is unity are related to regular continued fractions.

b) The values in the right-hand columns have a
simple relationship to the values in the left-hand columns;
can the reader see the correspondence before reading any further?
We have the identity
$$$1 - \bslash x↓1, x↓2, \ldotss , x↓n\bslash = \bslash 1,
x↓1 - 1, x↓2, \ldotss , x↓n\bslash ;\eqno (45)$
|qctr \6skip see exercise 9.
|qleft \qquad c) There is symmetry between left and right in
the first two columns: If $\bslash A↓1, A↓2, \ldotss , A↓t\bslash$
occurs, so does $\bslash A↓t, \ldotss , A↓2, A↓1\bslash $.
This will always be the case (see exercise 26).

d) If we examine all of the quotients in the table,
we find that there are 96 in all, of which ${39\over 96}$ =
40.6 percent are equal to 1, ${21\over 96}$ = 21.9 percent are
equal to 2, ${8\over 96}$ = 8.3 percent are equal to 3; this
agrees reasonably well with the probabilities listed above.

\subsectionbegin{The number of division steps} Let
us now return to our original problem and investigate $T↓n$,
the average number of division steps when $v = n$. $\biglp$See
Eq.\ (19).$\bigrp$\xskip Here are some sample values of
$T↓n:$
$$
|qctr \hfill n ⊗= 95 \xskip 96 \xskip 97 \xskip 98 \xskip 99
100 101 102 103 104 105\cr
\3skip \hfill T$↓n ⊗= 5.0 4.4 5.3 4.8 4.7 \xskip 4.6 \xskip
5.3 \xskip 4.6 \xskip 5.3 \xskip 4.7 \xskip 4.6\cr
\hfill n ⊗= 996 997 998 999 1000 1001 \cdots 9999 10000 10001\cr
\3skip \hfill T↓n ⊗= 6.5 \xskip 7.3 \xskip 7.0 \xskip 6.8 \xskip
6.4 \xskip 6.7 \cdots \xskip 8.6 \quad 8.3 \quad 9.1\cr
\3skip \hfill n ⊗= 49999 50000 50001 \cdots 99999 100000 100001\cr
\3skip \hfill T↓n ⊗= \xskip 10.6 \quad 9.7 \xskip 10.0 \cdots
\xskip 10.7 \quad 10.3 \quad 11.0\cr
\6skip$ Note the somewhat erratic behavior; $T↓n$ tends to be
higher than its neighbors when $n$ is prime, and it is correspondingly
lower when $n$ has many divisors. (In this list, 97, 101, 103,
997, and 49999 are primes; 10001 = 73 \cdot 137, 50001 = 3 \cdot
7 \cdot 2381, 99999 = 3 \cdot 3 \cdot 41 \cdot 271, 100001 =
11 \cdot 9091.) It is not difficult to understand why this happens;
if gcd($u, v) = d$, Euclid's algorithm applied to $u$ and $v$
behaves essentially the same as if it were applied to $u/d$
and $v/d$. Therefore, when $v = n$ has several divisors, there
are many choices of $u$ for which $n$ behaves as if it were
smaller.

Accordingly let us consider {\sl another} quantity,
$\tau ↓n$, which is the average number of division steps when
$v = n$ and when $u$ is {\sl relatively prime} to $n$. Thus
$$$\tau ↓n = {1\over \varphi (n)} \sum ↓{0≤m<n,??$gcd$(m,n)=1}
T(m, n).\eqno (46)$
|qctr \6skip It follows that
$$$T↓n = {1\over n} \sum ↓{d\rslash n} \varphi (d)\tau ↓d.\eqno
(47)$
|qctr \6skip Here is a table of $\tau ↓n$ for the same values
of $n$ considered above:
$$
|qctr \hfill n ⊗= 95 \xskip 96 \xskip 97 \xskip 98 \xskip 99
100 101 102 103 104 105\cr
\3skip \hfill \tau $↓n ⊗= 5.4 5.3 5.3 5.6 5.2 \xskip 5.2 \xskip
5.4 \xskip 5.3 \xskip 5.4 \xskip 5.3 \xskip 5. \cdot \cdot 9999
10000 10001\cr
\3skip \hfill \tau ↓n ⊗= 7.2 \xskip 7.3 \xskip 7.3 \xskip 7.3
7.3 7.4 \cdots 9.21 9.21 9.22\cr
\3skip \hfill n ⊗= 49999 50000 50001 \cdots 99999 100000 100001\cr
\hfill \tau ↓n ⊗= 10.58 \xskip 10.57 \xskip 10.59 \cdots 11.170 \xskip
11.172 \xskip 11.172\cr
\6skip$ Clearly $\tau ↓n$ is much more well-behaved than $T↓n$,
and it should be more susceptible to analysis. Inspection of
a table of $\tau ↓n$ for small $n$ reveals some curious anomalies;
for example, $\tau ↓{50} = \tau ↓{100}$ and $\tau ↓{60} = \tau
↓{120}$. But as $n$ grows, the values of $\tau ↓n$ behave quite
regularly indeed, as the above table indicates, and they show
no significant relation to the factorization properties of $n$.
If the reader will plot the values of $\tau ↓n$ versus ln $n$
on graph paper, for the values of $\tau ↓n$ given above, he
will see that the values lie very nearly on a straight line,
and that the equation
$$$\tau ↓n \approx 0.843$ ln $n + 1.47\eqno (48)$
|qctr \6skip gives a very good approximation to $\tau ↓n$.

We can account for this behavior if we return to
the considerations of Theorem W\null. Note that in Euclid's algorithm
as expressed in (15) we have
$$${V↓0\over U↓0} {V↓1\over U↓1} \ldotsm {V↓{t-1}\over U↓{t-1}}
= {V↓{t-1}\over U↓0} $,
|qctr \6skip since $U↓{k+1} = V↓k$; therefore, if $U = U↓0$
and $V = V↓0$ are relatively prime, and if there are $t$ division
steps, we have
$$$X↓0X↓1 \ldotsm X↓{t-1} = 1/U$.
|qctr \6skip Setting $U = N$ and $V = m < N$, we find that
$$ln $X↓0 +$ ln $X↓1 +\cdot \cdot \cdot +$ ln $X↓{t-1} = -ln
$N.\eqno (49)$
|qctr \6skip We know the approximate distribution of $X↓0, X↓1,
X↓2, \ldotss $, so we can use this equation to estimate
$$$t = T(N, m) = T(m, N) - 1$.
|qctr \quad Returning to the formulas preceding Theorem W\null, we
find that the average value of ln $X↓n$, when $X↓0$ is a real
number uniformly distributed in [0, 1), is
$$\6skip $\6skip \int ↑{1}↓{0}$ ln $xF↑\prime↓{n}(x)dx =
\int ↑{1}↓{0}$ ln $xf↓n(x)dx/(1 + x),\eqno (50)$
|qctr \12skip where $f↓n(x)$ is defined in (31). Now
$$$f↓n(x) = {1\over$ ln 2} + O(2↑{-n}),\eqno (51)$
|qctr \6skip using the facts derived earlier (see exercise 23);
hence the average value of ln $X↓n$ is, to a very good approximation,
|qleft β?????????|newtype 58320---Compute